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DARPA  SCIENTIFIC  REPORT  91-93 


During  these  three  years  of  activity  we  have  concentred  our  efforts  on  the 
following  subjects  : 


I  -  Regional  waves  propagation  and  attenuation  which  are  one  of  our  basic 
topics  for  several  years, 


II  -  Modeling  of  local  and  teleseismic  waves  generated  by  underground 
nuclear  explosions  in  the  NTS  and  in  the  Hoggar  Massif,  Sahara. 


Ill  -  First  data  processing  associated  with  a  mini-array  recently  built  in  the 
Center  of  France. 
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I  -  Regional  waves  propagation  : 


I  -  1  Crustal  wave  propagation  anomaly  across  the  Pyrenean  Range  : 

Following  our  previous  studies,  we  are  using  Lg  records  analysis  and  numerical 
modelling  of  Lg  propagation  to  find  out  to  what  extent  this  phase  can  be  seen  as  a 
marker  of  unidentified  structural  anomalies  within  the  crust. 

This  study  is  based  on  Lg  propagation  through  the  Pyrenean  range  from  earthquakes 
located  in  Spain. 

We  have  first  evaluated  the  mean  value  of  the  S-wave  quality  factor  for  Central-Spain. 
We  have  computed  simultaneously  the  seismic  station  responses  and  the  source 
functions.  The  correction  for  propagation  effects,  assuming  an  homogenous  attenuation 
and  the  theoretical  evaluation  of  the  Lg  excitation  lead  to  the  seismic  moment  of  each 
event.  The  moment  magnitude  we  obtained  fits  the  magnitude  proposed  by  the  local 
networks. 

This  gives  the  confirmation  of  the  Q  model  in  the  low  frequency  range  (l-5Hz).  As  we 
intended  to  compare  traces  of  different  Spanish  earthquakes  recorded  in  France  at 
different  epicentral  distances,  we  had  to  make  amplitudes  independent  of  propagation 
and  sources  effects.  Therefore  we  corrected  the  spectral  amplitudes  for  geometrical 
spreading,  anelastic  attenuation  and  normalized  them  to  equal  seismic  moment. 

We  then  plotted  the  records  as  a  function  of  group  velocity,  in  order  to  make  up  a  fan 
profile  along  the  Pyrenean  axis.  The  resulting  section  reveals  that  in  the  central  and  the 
eastern  parts  of  the  range,  neither  the  North-Pyrenean-Fault,  nor  the  Moho  jump 
deduced  from  seismic-refraction  experiments  and  vertical  seismics,  seem  to  affect  the 
Lg  propagation.  However,  there  is  an  extinction  of  the  Lg  phase  in  the  western  part  of 
the  chain.  The  lateral  extent  of  this  area  is  correlated  with  a  zone  of  positive  gravity 
anomaly,  probably  linked  to  the  presence  of  dense  material  of  mantle  origin. 

A  numerical  simulation  in  the  low  frequency  band  indicates  that  the  Moho  topography 
inferred  from  deep  seismic  soundings  does  not  explain  the  strength  of  the  observed 
attenuation.  Ray-tracing  seismograms  show  that,  at  high  frequency,  the  conclusion  is 
the  same. 
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The  attenuation  effect  of  the  structure  lateral  variation  should  not  be  so  strong.  We 
therefore  think  that  attenuation  of  guided  waves  is  not  basically  due  to  large  scale 
geometry  effects,  but  more  probably  to  local  properties  of  the  crustal  materials, 
possibly  apparent  attenuation  due  to  scattering  on  small  scale  heterogeneities. 

This  work  has  been  submitted  to  the  Geophysical  Journal  and  revised  in  march  1993. 


1-2  Calculation  of  synthetic  seismograms  in  a  laterally  varying  medium  by  the 
boundary  element  discrete  wavenumber  method  : 

We  present  a  theoretical  investigation  of  the  effect  of  lateral  crustal  heterogeneities  on 
the  propagation  of  seismic  waves.  We  study  particularly  how  Lg  waves  are  affected  by 
the  crossing  of  complex  crustal  structures.  The  work  is  carried  out  by  numerical 
simulation  and  the  method  of  calculation  is  based  on  boundary  integral  /  boundary 
element  techniques.  The  Green's  functions  necessary  to  the  implementation  of  these 
techniques  are  evaluated  in  the  frequency  -  wavenumber  domain  using  the  discrete 
wavenumber  method  and  reflectivity/transmissivity  matrices.  This  formulation  has  the 
advantage  over  purely  numerical  methods  of  allowing  the  calculation  of  the  wavefields 
over  large  distances  (several  hundred  kilometers)  corresponding  to  several  hundred 
times  the  wavelengths. 


We  investigate  how  crustal  waves  are  affected  by  the  crossing  of  faults  and  by  the 
presence  of  lateral  variations  in  crustal  thickness.  We  find  that  the  Lg  waveforms  are 
extremely  sensitive  to  the  presence  of  lateral  heterogeneities  along  their  path,  but  that 
their  energy  stays  remarkably  stable  as  long  as  the  heterogeneities  encountered  are  not 
too  strong.  On  the  basis  of  the  theoretical  results  obtained,  the  observations  of  Lg  wave 
extinction  or  strong  attenuation  over  specific  continental  paths  are  difficult  to  explain 
by  purely  geometric  diffraction  and  scattering. 

The  numerical  simulation  results  also  show  that  the  back  scattered  Lg  wavefield 
induced  by  lateral  variations  in  crustal  thickness  can  be  strong.  This  suggests  that  in 
real  observations  part  of  the  Lg  wave  coda  is  due  to  back  scattering  and  that  the  level  of 
coda  present  might  be  related  to  the  roughness  of  the  Moho. 


3 
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SUMMARY 

Lg  records  analysis  and  numerical  modelling  of  Lg  propagation  are  used  to  find  out  to  what  extent  this  phase 
can  be  seen  as  a  marker  of  unidentified  structural  anomalies  in  the  crust  This  study  is  based  on  Lg 
propagation  through  the  Pyrenean  tange  from  earthquakes  located  in  Spain. 

We  have  first  evaluated  the  mean  value  of  the  5-wave  quality  factor  for  Central  Spain  We  have  computed 
simultaneously  the  seismic  station  responses  and  the  source  functions  The  correction  for  propagation  effects, 
assuming  a  homogeneous  attenuation  and  the  theoretical  calculation  of  the  Lg  excitation,  lead  to  the  seismic 
moment  of  each  event.  The  moment  magnitude  obtained  correlates  well  with  the  magnitude  proposed  by  the 
local  networks.  This  gives  a  confirmation  of  the  0  model  in  the  low  frequency  range  (l-5Hz).  As  we  intended 
to  compare  traces  of  different  Spanish  earthquakes  recorded  in  France  at  different  epicentral  distances,  we  had 
to  make  amplitudes  independent  of  propagation  and  source  effects.  Therefore,  we  corrected  the  spectral 
amplitudes  for  geometrical  spreading,  anelastic  attenuation  and  normalized  them  to  equal  seismic  moment. 

We  then  plotted  the  records  as  a  function  of  group  velocity,  in  order  to  make  up  a  fen  profile  along  the 
Pyrenean  axis  The  resulting  section  reveals  that  in  the  central  and  the  eastern  parts  of  the  range,  neither  the 
North-Pyrenean-Fauh,  nor  the  Moho  jump  deduced  from  seismic-refraction  experiments  and  vertical  seismics, 
seem  to  affect  Lg  propagation.  However,  there  is  an  extinction  of  the  Lg  phase  in  the  western  part  of  the  chain. 
The  lateral  extent  of  this  area  is  correlated  with  a  zone  of  positive  gravity  anomaly,  probably  linked  to  the 
presence  of  dense  material  of  mantle  origin.  A  numerical  simulation  in  the  low  frequency  band  indicates  that 
the  Moho  topography  inferred  from  deep  seismic  soundings  does  not  explain  the  strength  of  the  observed 
attenuation.  Ray-tracing  seismograms  show  that,  at  high  frequency  the  conclusion  is  the  same.  The  attenuation 
effect  of  the  structure  lateral  variation  should  not  be  so  strong.  We,  therefore,  think  that  attenuation  of  guided 
waves  is  not  due  to  large-scale  geometry  effects,  but  is  due  to  local  properties  of  the  crustal  materials,  possibly 
apparent  attenuation  due  to  scattering  cm  small-scale  heterogeneities. 

Key  words  :  Lg  waves,  Pyrenees,  quality  factor,  synthetic  seismograms. 

INTRODUCTION 

In  the  range  between  150  and  several  thousand  kilometres,  crustal  waves  are  dominant  on  short-period 
seismograms  in  continental  areas.  Thus,  the  Lg  phase,  which  consists  of  5  waves  trapped  in  the  crust,  is  the 
major  phase  observed  on  regional  records. 

Lg  amplitude  is  known  to  be  sensitive  to  important  changes  in  crustal  structure:  propagation  paths  through 
oceanic  crust  are  the  origin  of  high  attenuation  or  extinction  of  the  Lg  phase,  as  found  in  the  early  analysis  of 
this  phase  (Press  &  Ewing  1952;  Bath  1954).  Zones  of  strong  local  weakening  of  Lg  also  exist  in  continental 
domains.  Such  observations  have  been  reported  in  the  Himalayan  Belt  (Ruzaikin  et  al.  1977),  in  the  North  Sea 
Graben  (Kennett  et  al.  1985),  in  the  Eurasian  Shield  (Baumgardt  1991)  and  in  the  south-western  part  of  the 
Alpine  range  (Campillo  et  al.  1993). 

ft  seems  that  some  geological  features  partially  or  completely  stop  crustal  guided  wave  propagation  Our 
purpose  is  to  investigate  what  structure  may  produce  these  Lg  amplitude  variations  and  blockage  effects.  Is  it  a 
large-scale  geometry  result  of  a  local  attenuation  or  scattering  effect?  As  there  already  exists  evidence  of  an 
Lg  propagation  anomaly  in  the  Western  Alps,  our  study  has  been  aimed  at  the  analysis  of  Lg  propagation 
across  the  Pyrenean  range,  in  order  to  improve  our  knowledge  about  the  influence  of  such  orogenic  structures. 
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DATA  SPECIFICATIONS 


For  this  study,  we  have  used  Spanish  and  French  data  obtained  by  the  IGN  (Instituto  Geografico  Nacional) 
array  and  the  LDG  (  Laboratoire  de  Detection  Geophysique)  network  respectively  All  seismic  stations  are 
short -period  vertical  instruments  with  natural  period  of  Is.  The  characteristics  of  the  networks  are  described  m 
an  IGN  report  (IGN  1991)  and  by  Nicolas  et  al.  (1982)  for  the  LDG  network.  We  have  examined  records 
available  from  these  two  networks  and  we  confirmed  that  Lg  waves  effectively  propagate  in  France  and  in 
Spain,  on  both  sides  of  the  Pyrenean  range.  As  we  intended  to  study  die  propagation  from  natural  sources 
located  in  Spain,  to  see  whether  or  not  Lg  waves  are  able  to  cross  the  Pyrenees,  we  had  to  determine  the 
anelastic  attenuation  of  the  crustal  phases  in  Spain 

The  data  set  consists  of  12  earthquakes  with  hypocentres  in  the  crust,  recorded  in  Spain.  Their  parameters  are 
described  in  Table  1 .  As  is  the  case  in  most  continental  areas.  Lg  is  the  dominant  phase  on  most  of  the  records. 
The  source-receiver  pairs  have  been  chosen  in  such  a  way  that  the  propagation  paths  were  sampling  the  central 
part  of  Spain,  considered  to  be  a  stable  continental  area.  We  did  not  take  into  account  stations  located  in 
southern  Spain  because  we  observed  strong  attenuation  in  the  region  of  Gaudalquivir  sedimentary  basin  Fig  1 
shows  the  paths  used  in  this  part  of  the  study 

EVALUATION  OF  THE  5-WAVE  MEAN  QUALITY  FACTOR 

To  evaluate  the  5-wave  mean  quality  factor,  we  used  Lg  phases  that  consist  of  multiply  reflected  5  waves 
Calculation  of  the  crustal  quality  factor  O  is  done  by  computing  the  spectral  density  per  time  unit  as  a  function 
of  epicentral  distance,  in  the  time  window  providing  the  largest  amplitudes,  i.e.  in  the  group  velocity  window 
3.6-3 .2  km  s'1 


For  each  station  /  and  each  earthquake  j,  this  amplitude  can  be  written  in  the  form  : 

Aij  (f  d)  =  Sj( f)  .  AA( l  d)  .  E(d)  .  STi(f),  (1) 


f  representing  the  frequency  and  d  the  epicentrral  distance.  Sj(i)  is  the  source  contribution,  proportional  to  the 
seismic  moment  at  low  frequency  AA(f,  d)  is  given  by 


exp 


-/rfd 

yOVm  J 


(2) 


and  represents  the  anelastic  attenuation  for  a  phase  propagating  with  a  mean  velocity  Q\’m.  E(d)  denotes  the 
geometrical  spreading  in  the  time  domain  needed  to  correct  the  spectral  density  per  time  unit.  It  is  taken  in  the 
form: 


E(d)  =  d-083  (Campillo,  Bouchon  &  Massinon  1984)  (3) 

ST/(f)  is  the  station  response,  corresponding  to  the  effects  of  local  geology  structure. 


This  equation  is  solved  by  an  iterative  process  given  in  detail  by  Campillo,  Plantet  &  Bouchon  (1985).  Only 
the  data  obtained  for  earthquakes  that  were  recorded  by  at  least  four  stations,  were  used.  We  first  took  ST  =  1 
and  evaluated  the  mean  value  of  O  by  least-square  fitting  for  each  earthquake.  We  then  computed  the  mean 
value  of  Q  for  the  entire  set  of  events.  ST  could  be  evaluated  by  computing  a  simple  residue  at  each  station. 
This  process  converged  after  a  few  iterations.  The  root-mean-square  residue  between  the  observed  amplitudes 
and  the  ones  predicted  by  eq.  1  was  computed  to  check  the  convergence  and  the  stability.  The  shape  of  the 
displacement  spectrum  was  obtained  after  deconvolution  of  the  instumental  response  of  the  IGN  network 
stations.  Thus,  we  got  the  sources-displacement  spectra  inms. 

We  found  0(f)  in  the  form  : 

O  =(330±30)/°",±oos  (4) 


This  result  is  close  to  the  mean  crustal  quality  factor  computed  for  central  France  (Campillo  et  al. 
1985) 
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ESTIMATION  OF  THE  SEISMIC  MOMENTS 


We  computed  the  seismic  moment  using  the  theoretical  excitation  of  Lg  for  a  point  source  dislocation.  We 
neglected  the  radiation  pattern  since  Lg  is  a  superposition  of  S  waves  leaving  the  source  in  a  wide  range  of 
take-off  angles.  To  perform  the  estimation,  we  computed  synthetic  seismograms  in  a  flat  layered  medium 
corresponding  to  the  crustal  structure  of  central  Spain.  Table  2  summarizes  the  characteristics  of  the  model. 

The  theoretical  calculation  was  performed,  for  a  seismic  moment  of  1,  using  the  discrete  wavenumber 
representation  (Bouchon  19Sl),  combined  with  the  Kennett  propagation  technique  (1983).  We  evaluated  the 
Lg  spectral  density  from  the  synthetics  in  exactly  the  same  way  as  for  the  data.  We  can,  therefore,  obtain  the 
theoretical  value  of  S  of  eq.  1  for  a  unit  moment.  For  frequencies  lower  than  the  comer  frequency,  we  can 
relate  the  seismic  moment  to  the  observed  amplitudes,  corrected  for  spreading  and  attenuation.  To  test  the 
accuracy  of  our  results  we  have  calculated  the  Lg  magnitudes  from  each  seismic  moment  value,  from 
which  we  deduced  Mw  as  defined  by  Kanamori  (1977)  is  given  by: 

logA/o  =  \.5Mw  +  9  {Mo  in  Nm).  (4) 

Table  1  presents  the  local  magnitudes  reported  in  the  bulletins  and  the  values  of  our  magnitudes.  We  can 
see  that  the  seismic  moments  measured  from  the  Lg  phase  vary  consistently  with  the  magnitudes  proposed  by 
the  local  networks.  Fig  2  shows  the  source  spectra  which  allow  us  to  find  out  the  seismic  moments  and  the 
magnitude  values.  We  plotted  all  spectra  on  a  log-log  diagram,  in  order  to  compare  their  shape  and  their 
dimension.  One  can  see  that  the  comer  frequency  measured  chi  our  spectra  varies  between  3  an  8  Hz. 
Considering  the  events  with  seismic  moment  between  1013  and  1016  Nm  i.e.  Mw  <5,  and  assuming  a  self¬ 
similar  model,  fc  should  be  in  the  range  0.6-7  Hz  for  a  100  to  200  bars  stress  drop.  Beyound  fc,  the  observed 
high-frequency  spectra  decrease  as  w  "2. 

The  Lg  magnitudes  we  found  are  systematically  smaller  than  those  given  by  the  local  networks  from  direct 
waves.  We  have  tested  that  the  crustal  model  used  for  the  numerical  simulation  of  Lg  propagation  has  no 
significant  influence  on  the  seismic  moment  values.  On  the  other  hand,  changes  in  the  parameters  of  fault 
geometry  cause  variations  in  the  Lg  magnitudes.  However,  azimuthal  dependence  is  weaker  for  Lg  than  for 
direct  waves  (Campillo  1990)  because  this  phase  is  made  up  of  a  superposition  of  rays,  sampling  a  wide  range 
of  take-off  angles.  The  good  agreement  of  the  linear  correlation  we  obtained  between  Lg  moment  magnitudes 
and  local  magnitudes  confirms  the  possibility  of  using  the  source  spectra  to  correct  the  amplitudes. 

PROPAGATION  ANALYSIS  THROUGH  THE  PYRENEES 

Most  of  the  earthquakes  whose  records  were  used  for  the  Q  calculation  did  not  provide  data  of  adequate 
quality  at  the  French  stations  because  of  the  great  epicentral  distances.  Only  four  events  located  in  northern 
Spain  provide  a  large  enough  signal-to-noise  ratio:  Sotos,  Cucalon  and  the  two  closely  spaced  events  Camera 
and  Amedo  (see  Fig.  1). 

First  we  looked  at  the  records  obtained  at  the  French  station  EPF,  which  is  located  in  the  Pyrenees,  north  of 
the  North  Pyrenean  Fault  (NPF,  Fig.  3).  This  zone  of  deep  subvertical  faults  is  a  major  structural  feature  of 
the  mountain  range,  clearly  apparent  in  the  oriental  and  in  the  central  Pyrenees.  In  the  western  Pyrenees,  it  is 
assumed  that  the  discontinuity  is  overlain  by  Cretaceous  sediments.  The  NPF  constitutes  the  limit  of  the  North 
Pyrenean  Zone  and  the  axial  zone  of  the  range.  The  available  information  concerning  the  tectonics  of  the 
region  shows  that,  in  the  central  part  of  range  (Fig.  4),  an  underthrust  of  the  Iberic  crust  to  the  North  takes 
place  along  the  NPF  (Pinet  et  al.  1987;  Choukroune  et  al.  1989).  On  the  other  hand,  the  oceanic  lithosphere  in 
the  Gulf  of  Biscay  dives  under  the  Iberic  plate  (Boillot  et  al.  1971). 

For  the  earthquakes  we  examined,  the  seismograms  recorded  at  EPF  have  the  same  typical  shape  of 
continental  short-period  records  as  observed  at  all  nearby  Spanish  stations:  the  crustal  wave  Pg  and  Lg  are 
dominant.  This  is  illustrated  in  Fig.  5  (a  and  b):  the  Lg  phase  produces  a  clear  onset  at  3.5  km  s'1  .  However, 
the  records  obtained  at  stations  in  central  France  show  very  different  characteristics.  The  corresponding  paths 
are  plotted  in  Fig.  6.  In  the  case  where  the  path  crosses  foe  western  part  of  foe  Pyrenean  Chain,  foe  larger 
arrivals  on  foe  seismograms  have  a  group  velocity  higher  than  Lg  (about  4.2  km  s'*)  and  consist  of  foe  mantle 
wave  Sn  (Fig.  5c).  The  Lg  phase  vanishes  along  this  path.  For  a  path  that  crosses  the  central  part  of  foe  chain, 
Sn  and  Lg  exhibit  similar  amplitudes  (Fig.  5d).  These  seismograms  suggest  that  foe  regional  phases  are 
affected  very  differently  when  they  cross  foe  Pyrenees  at  different  locations  along  foe  axis  of  foe  range. 
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In  order  to  verify  this  effect  we  plotted  a  senes  of  seismograms  as  a  function  of  the  location  where  the  waves 
cross  the  Pyrenees.  We  measured  the  locations  of  the  crossing  along  an  axis  from  Bilbao  to  Perpignan  shown 
as  a  thick  solid  line  in  Fig.  6.  The  identification  of  the  various  regional  phases  in  the  data  set  is  simplified  by 
plotting  the  records  as  a  function  of  group  velocity.  We  consider  the  records  in  France  from  a  series  of 
earthquakes  in  northern  Spain  (Table  3  and  Fig.  6).  The  records  are  bandpassed  between  1  and  5  Hz. 
Amplitudes  are  corrected  for  geometrical  spreading,  anelastic  attenuation  and  normalized  to  equal  seismic 
moment.  The  effect  of  attenuation  is  crudely  removed  from  the  whole  seismogram  by  correcting  the  spectra 
with  a  filter  defined  from  our  results  of  the  quality  factor  for  S  waves  (eq.  4).  The  seismic  moment  used  for  the 
normalization  is  obtained  from  the  least-square  regression  of  the  spectral  amplitude  of  Lg  with  distance  for  the 
LDG  stations.  In  case  of  a  strong  effect  of  attenuation  of  Lg  at  the  crossing  of  die  Pyrenees,  the  moment  is 
underestimated.  This  results  in  an  unrealistically  high  amplitude  of  the  other  seismic  phases  such  as  Sn. 

Figure  7  presents  the  section  obtained  from  diese  Spanish  earthquakes  recorded  in  France.  One  can  see  that  the 
sampling  is  not  homogeneous  because  of  the  poor  distribution  of  earthquakes  in  northern  Spain  correctly 
recorded  in  France.  However,  we  notice  that  the  waveform  is  very  different  in  the  east  and  in  the  west  side  of 
the  Pyrenean  range.  In  the  east,  the  maximum  amplitudes  appear  for  group  velocities  between  3.5  an  3.0 
km  s"1  these  waves  are  clearly  Lg  phases.  In  the  west,  the  largest  amplitudes  are  seen  between  4.5  and  4  0  km 
s'1  .  These  phases  are  Sn  waves  corresponding  to  upper  mantle  propagation. 

As  has  already  been  noticed,  the  way  we  perform  the  normalization  of  seismic  moment  can  explain  why  the 
amplitudes  of  Pn  and  Sn  waves  in  the  western  part  appear  much  higher  than  those  in  die  eastern  part.  Another 
reason  is  the  way  we  corrected  amplitudes  with  distance.  The  amplitude  corrections  for  spatial  attenuation  are 
strictly  valid  only  in  the  group  velocity  window  3.2-3. 6  km  s*l,  i.e.  for  die  Lg  wave.  Phases  for  which  group 
velocity  is  greater  than  the  Lg  wave  velocity  are  all  the  more  amplified  since  epicentral  distance  and  frequency 
increase.  Even  if  the  amplitude  values  are  not  exact  along  die  entire  record,  this  figure  illustrates  die  Lg 
blockage  and  the  relative  Sw-wave  amplification  in  the  western  Pyrenees.  The  Sn  waves  cannot  be  observed  in 
the  eastern  Pyrenees  because  their  amplitudes  are  much  smaller  than  the  ones  of  Pg  and  Lg.  If  we  accept  that 
Sn  is  not  sensitive  to  crustal  structure  and  that  the  mantle  is  almost  homogeneous  beneath  this  region,  the 
variation  of  the  ratio  of  amplitude  Sn  to  Lg  at  a  given  distance  is  a  crude  measure  of  the  variation  of 
attenuation  in  the  crust.  Fig.  7  indicates  that  between  eastern  and  western  parts  of  the  Pyrenees,  the  attenuation 
effect  changes  by  more  than  one  order  of  magnitude. 

However,  for  all  these  earthquakes,  the  Lg  phase  can  be  observed  at  EPF  station.  The  presence  of  the  Lg  phase 
in  the  east  and  at  EPF  indicates  that  the  Lg  blockage  is  not  simply  associated  with  the  NPF,  which  lies  along 
the  entire  eastern  and  central  part  of  the  range.  On  the  other  hand,  the  Lg  phase  vanishes  when  it  crosses  the 
western  part  of  the  mountain  chain  and  the  energy  seems  to  propagate  mostly  as  the  Sn  wave.  This  anomaly 
seems  to  have  a  lateral  extent  of  about  100  km.  It  is  noteworthy  that  this  zone  of  attenuation  corresponds  to  a 
zone  of  strong  positive  gravity  anomaly  called  the  'Labourd  anomaly1,  whose  origin  is  not  definitively  known. 
The  areas  with  positive  Bouguer  anomaly  are  shown  in  grey  in  Fig.  6.  A  similar  anomaly  of  propagation  of 
the  crustal  phase  Lg  has  been  reported  in  the  western  Alps  (Campillo  et  al.  1993),  correlated  with  a  positive 
Bouguer  anomaly.  Crustal  materials  of  deep  origin  might  be  associated  with  both  features. 

NUMERICAL  SIMULATIONS  IN  MULTILAYERED  MEDIA  WITH  IRREGULAR  INTERFACES 

We  attempted  to  understand  why  Lg  is  not  observed  through  the  western  part  of  the  range.  To  investigate  the 
influence  of  crustal  geometry  on  Lg  propagation,  we  performed  some  numerical  modelling  in  laterally 
heterogeneous  models  for  the  SH  case.  Since  Lg  consists  of  a  superposition  of  post-critically  reflected  S  waves, 
the  SH  case  must  present  most  of  the  effect  of  the  large-scale  lateral  variations  of  Moho  depth.  The  calculation 
method  combines  the  discrete  wavenumber  Green's  function  representation  with  boundary-integral  equation 
techniques  (Campillo  &  Bouchon  1985;  Bouchon,  Campillo  &  Gaffet  1989).  The  wavefield  produced  by  the 
interfeces  is  considered  to  be  equivalent  to  the  radiation  of  body  forces  distributed  along  the  boundaries.  The 
inversion  of  a  propagator  matrix  was  performed  for  each  interface,  so  that  the  computation  time  and  the 
memory  required  varies  only  linearly  with  the  number  of  interfeces.  As  many  seismic  experiments  revealed  the 
presence  of  a  10  km  Moho  jump  between  the  North  Pyrenean  Zone  and  the  axial  zone  (Him  et  al.  1980;  Roure 
et  al.  1989),  we  have  first  designed  a  model  with  simple  change  in  crustal  thickness  (Fig  8a),  the  Moho  being 
deeper  on  the  Spanish  side  of  the  Pyrenees.  There  is  no  attenuation  in  this  model:  only  the  topography  of  the 
Moho  is  taken  into  account.  7 


The  synthetic  seismograms,  with  a  Ricker  wavelet  of  1.5  s  period  as  the  source  function,  are  shown  in 
Fig.  8(b).  The  maximum  frequency  reached  is  \  Hz.  The  reduction  velocity  is  3.5  km  s'1,  which  is  die  5-wave 
velocity  chosen  for  the  crust.  One  can  see  the  successive  reflection  branches,  which  constitute  the  reflected 
energy  forming  the  Lg  phase.  The  Sn  head-wave  branch  appears  for  distances  greater  than  250  km  with  small 
amplitudes. 

We  can  check  that,  at  low  frequencies,  the  Moho  jump  does  not  produce  a  notable  effect  on  wavesforms  and 
amplitudes.  One  can  notice  a  weak  perturbation  above  the  Moho  jump,  but  amplitudes  are  as  large  beyond  the 
jump  as  they  are  ahead  of  it.  The  use  of  boundary  integral  equations  is  limited  to  relatively  low  frequencies 
simply  because  of  the  high  computation  time  required  by  this  quasi-exact  approach.  In  order  to  check  the 
validity  of  our  conclusions  at  higher  frequency,  we  used  the  ray  theory  to  compute  an  ’infinite  frequency’ 
response.  The  calculations  were  made  with  the  same  model  of  Moho  topography.  Details  of  the  computations 
with  the  paraxial  ray  approach  and  of  the  comparison  with  the  boundary-integration  equation  method  results 
are  given  in  Gibson  &  Campillo  (1993).  The  synthetics  obtained  (Fig.  8c)  are  very  similar  to  those  of 
boundary  integral  equations  and  confirm  that  the  Moho  step  will  only  have  weak  effect  on  the  guided  wave 
amplitude,  whatever  the  frequency  band  is. 

From  these  numerical  tests  we  concluded  that  the  Moho  jump  cannot  be  the  reason  for  the  vanishing  Lg  wave. 
This  is  strongly  supported  by  the  observation  that  Lg  waves  are  not  attenuated  in  the  Eastern  Pyrenees.  We 
must,  therefore,  examine  in  more  detail  the  influence  of  the  particular  structure  of  the  western  part  of  the 
chain. 

Figure  9  (a)  describes  the  second  model,  which  includes  results  of  recent  deep  seismic  investigations  conducted 
in  the  western  Pyrenees  (preliminary  interpretation  of  the  structure  of  the  lithosphere  along  the  Pyrenees- 
Arzacq  Ecors  profile,  M.  Daignieres,  private  communication).  These  experiments  suggest  a  zone  of 
anomalously  high  velocity  in  the  crust,  in  addition  to  the  Moho  jump.  We,  therefore,  investigated  the  influence 
of  such  a  structure. 

The  synthetic  seismograms  from  this  second  model  ate  shown  in  Fig.  9(b).  One  can  notice  that  the  Sn  wave  is 
stronger  than  in  the  previous  model.  Therefore,  as  observed  on  the  real  data,  the  Sn/Lg  amplitude  ratio  is 
higher  for  the  rays  travelling  through  such  a  zone.  But  the  Lg  wave  is  not  much  affected:  the  amplitude  of  the 
Lg  wave  train  is  still  large  beyond  the  velocity  anomaly  area.  Therefore,  the  geometry  of  this  structure  can  not 
account  for  the  almost  complete  extinction  of  the  Lg  phase  observed  on  the  data. 

This  energy  blockage  must  then  be  explained  by  the  local  properties  of  the  crustal  material  rather  than  by  the 
large-scale  geometry  of  the  crust-mantle  structure.  The  theoretical  result  is  in  a  good  agreement  with  the 
observation  that  Lg  can  propagate  across  the  eastern  part  of  the  mountain  range  where  the  jump  of  the  Moho 
between  north  and  south  is  present  as  in  the  western  part.  The  high  velocities  observed  in  this  area  where  Lg 
disappears  may  be  due  to  the  presence  of  lower  crustal  blocks,  brought  up  to  the  surface  during  the 
compression  phases  of  the  orgeny,  or  due  to  the  presence  of  slices  of  mantle  materials.  Both  interpretations 
would  agree  with  existence  of  a  positive  Bouguer  anomaly.  As  the  lower  crust  is  known  to  be  layered  and  very 
reflective  in  many  parts  of  western  Europe  (see  Mooney  &  Brocher  1987,  for  a  review),  both  interpretations 
imply  an  increased  heterogeneity  of  the  crust  in  this  region  with  respect  to  neighbouring  areas.  An  enhanced 
scattering  of  seismic  waves  by  this  heterogeneity  may  be  the  case  of  the  strong  attenuation  of  crustal  phases 
observed  in  the  zone  where  the  gravity  anomaly  indicates  intrusion  of  material  of  deep  origin  into  the  upper 
crust. 


CONCLUSION 

We  combined  observations  and  numerical  simulations  in  order  to  study  crustal  wave  propagation  through  the 
Pyrenees.  We  first  computed  a  mean  crustal  value  of  the  5-wave  quality  factor  for  central  Spain  to  evaluate 
intrinsic  attenuation  along  the  paths.  O  has  been  found  in  the  form 

0  =  330f°-51 
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The  seismogram  analysis  shows  that  in  the  central  and  eastern  part  of  the  mountain  chain  the  North  Pyrenean 
Fault  does  not  block  Lg  propagation  and  that  the  Moho  jump  found  in  this  region  does  not  block  either.  In 
spite  of  the  feet  that  the  jump  of  the  Moho  is  present  all  along  the  mountain  range,  a  localized  zone  of 
attenuation  exists  in  the  western  part  of  the  Pyrenees,  correlated  with  a  positive  Bouguer  anomaly.  As  similar 
observations  were  made  in  the  Alps  in  the  Ivrea  region  and  as  numerical  modellings  show  that  geometrical 
effects  do  not  explain  the  observed  extinction  of  Lg  waves,  it  seems  that  the  general  conclusion  can  be  drawn 
that  strong  attenuation  of  guided  waves  is  probably  due  to  local  crustal  properties.  Scattering  by  small-scale 
heterogeneities,  such  as  lower  crust  or  mantle  slices,  may  be  fee  cause  of  strong  attenuation  in  fee  frequency 
range  considered  here.  This  interpretation  is  coherent  wife  fee  observation  of  high  seismic  velocity  and  fee 
position  of  a  Bouguer  anomaly  in  these  regions. 
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TABLE  1 


Event 

Date 

Latitude 

°N 

Longitude 

°W 

H 

Km 

Ml 

LDG 

MI 

IGN 

Mw-Lg 

Ayamonte 

20  12  89 

23.286 

-  07.328 

23. 

5.4 

5.0 

44.48 

Grandola 

23  09  89 

38.291 

-  08.655 

18. 

4.0 

38.40 

Cadix 

25  10  87 

36.361 

-  06.501 

66. 

4.8 

4.2 

36.40 

Alicante 

02  12  88 

38.200 

+  00.315 

20. 

3.8 

3.8 

33.37 

Sotos 

07  1 1  87 

43.000 

-03.812 

15. 

3.9 

3.8 

35.37 

Vilanova 

10  08  86 

41.128 

-  07.213 

12. 

4.0 

4.3 

36.39 

Laza 

30  08  89 

42.105 

-07.516 

13 

3.7 

3.9 

35.37 

Nazare 

31  03  89 

39.601 

-  09.493 

25? 

3.7 

3.5 

33.35 

Camera 

20  09  87 

42.138 

-  02.476 

05. 

3.5 

3.6 

34.35 

Aldea 

24  05  88 

38.400 

-  03.406 

12. 

3  8 

34.37 

Cucalon 

06  07  87 

40.951 

-01.001 

05. 

3.5 

3.3 

31.33 

Amedo 

29  11  88 

42.146 

-02.112 

05. 

3.1 

3.2 

29.31 

TABLE  2 

Crustal  model  used  for  the  numerical  calculations 

Layer  s 

Thickness 

(km) 

P.  Wave 

(km/sec) 

S.Wave 

Velocity 

(km/sec) 

Density 

(g/cm3) 

P.Wave 

Q 

S.wave 

Q 

2 

3.3  . 

2.5 

2.1 

10000 

10000 

5 

6.1 

3.48 

2.8 

10000 

10000 

4 

5.6 

3.18 

2.7 

10000 

10000 

13 

6.4 

3.58 

2.9 

10000 

10000 

7 

6.85 

3.90 

3.0 

10000 

10000 

8.05 

4.45 

3.3 

10000 

10000 

ll 


TABLE  3 


Date 

Time 

Latitude 

°N 

Longitude 

°W 

Ml 

LDG 

07  04  86 

23h32'07 

42.91 

-  1.90 

3.1 

27  10  86 

06h48'12 

39.88 

- 1.12 

3.5 

24  08  87 

18h43'06 

41.01 

1.49 

4.4 

20  02  88 

16h38'48 

42.40 

1.48 

3.8 

16  03  88 

21hl9'02 

42.38 

2.17 

3.8 

26  07  90 

16h29'32 

42.37 

-  1.29 

3.7 

20  08  90 

12h22'33 

40.26 

-0.92 

3.4 
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Figure  Captions 

Figure  1:  Locations  of  earthauake  epicenters  istarss  ano  IGN  networK  stations  (circles).  On  the  rr.aD  are  also 
reported  the  pains  usea  in  the  calculation  of  tre  mean  vaiue  of  crustal  quality  factor  in  Soam. 

F;gure  2:  Source  disDlacement  spectra  piottec  tor  eacn  earthouaKe  m  m.s. 

Figure  3:  Tectonic  map  of  the  Pyrenean  region.  7 he  station  EPF  ;s  reoortea  on  the  Figure.  NPF:  Norm  Pyrenean 
fault. 

Fioure  Cross  section  showing  the  crust  geometry  ar.c  the  Mono  tcpograpny  from  the  Ebro  oasin  in  the  South  to 
the  Aauitame  oasin  in  the  North  in  the  centra;  oar:  of  me  mountain  range  (from  Roure  et  ai..  1989). 

Figure  5:  Short  period  records  obtained  for  earthauakes  "Camero"  and  "Cucalon"  at  station  EPF  (a  anc  b)  and  in 
central  France  ic  and  d).  The  corresocncina  oatns  are  oiottec  in  Figure  6.  The  grouo  velocity  are  maicatea  in 
orcer  to  maKe  easier  the  icentificaticn  of  the  cifferent  phases. 

Figure  6:  MaD  showing  the  location  of  earthquakes  in  Spam-  (black  dots)  and  the  seismic  stations  in  France 
(circles)  used  to  study  regional  phases  crossing  the  Pyrenees.  The  lines  correspond  to  the  path  of  the 
seismograms  snown  in  Figure  5.  The  neavv  line  r.cicates  the  axis  usea  to  locate  the  crossing  of  the  cnain  for  the 
cifferent  pains.  The  grey  zones  mcicate  a  positive  Bouguer  gravity  anomaly. 

Figure  7:  Seismograms  obtained  in  France  for  earthquakes  in  Soain  piotted  as  function  of  group  velocity  to  allow  a 
direct  comparison  of  the  traces.  T'ne  horizontal  axis  reoresents  the  position  of  the  crossing  of  eacn  path  with  the 
line  Bilbao-Perpignan  shown  in  Figure  6.  Amputuces  are  corrected  from  propagation  effects  and  normalized  to 
equal  seismic  moment  using  the  propagation  parameters  obtained  in  central  Spain  and  central  France. 

Figure  8:  Influence  of  the  Mohc  topography  on  Lg  propagation,  a)  model  used  with  a  simple  variation  of  the  Mono 
depth,  b)  syntnetic  seismograms  obtained  using  the  boundary  integral  eouation  method.  c)synthetic  seismograms 
obtained  using  the  asymtotic  ray  theory.  Ail  syntnetics  are  piotteo  using  a  reduction  velocity  of  3.5  km/s. 

Figure  9:  Influence  of  the  Moho  tcpograohy  on  Lg  prcoagaticn.  a)  mooel  used  with  a  Moho  jump  and  introduction  of 
high  velocity  cccies  in  tne  crust,  b)  synthetic  seismograms  oDtamed  using  the  ocunoary  integral  equation  method. 
Synthetics  are  clotted  using  a  reauction  velocity  of  3.5  km/s. 
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Abstract 


We  investigate  the  propagation  of  Lg  waves  in  laterally-varying  crustal 
structures  by  numerical  simulation.  The  method  of  calculation  is  formulated  in 
terms  of  boundary  integral  equations  where  the  Green’s  functions  are  evaluated  by 
wavenumber  summation.  The  approach  is  well  suited  to  study  the  propagation  of 
seismic  waves  in  a  layered  medium  where  the  interfaces  are  fiat  in  some  regions  and 
irregular  in  other  regions.  We  investigate  the  effect  of  a  crustal  fault  with  vertical 
offset  and  study  the  case  of  a  lateral  change  in  crustal  thickness.  The  results  show 
that  the  Lg  wave  amplitude  is  only  slightly  affected  by  the  presence  of  these 
heterogeneities.  They  confirm  the  robustness  of  Lg  wave  propagation  in  presence  of 
lateral  heterogeneities  observed  in  other  numerical  simulations.  They  show  that 
large  scale  geometric  features  of  ihe  crust  cannot  account  alone  for  the  strong 
attenuation  of  Lg  waves  observed  in  many  regions.  The  results  also  suggest  a 
possible  relation  between  the  level  of  Lg  wave  coda  and  the  degree  of  roughness  of 
the  Moho.  They  further  indicate  the  importance  of  back  scattering  and  suggest  a 
possible  use  of  the  back  scattered  wive  field  to  map  strong  crustal  heterogeneities. 
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Introduction 


Over  the  last  two  decades,  numerical  simulation  techniques  have  become  an 
increasingly  important  tool  in  studying  the  propagation  of  seismic  waves  in  complex 
geological  structures.  Several  methods  of  elastodynamic  calculations  have  been 
developed  to  this  effect.  The  most  widely  used  techniques  may  be  classed  in  three 
groups:  ray  methods,  finite-diffenmce  finite-element  techniques,  and  boundary 
integral  equation  methods. 

These  numerical  simulation  studies  have  been  performed  for  different  types  of 
geological  structures.  Among  the  most  extensively  studied  geological  objects  are 
sedimentary  basins  and  alluvial  valleys  (Aki  and  Larner,  1970;  Trifunac,  1971; 
Bouchon  and  Aki,  1977a;  Hong  and  Helmberger,  1978;  Sanchez-Sesma  and  Esquivel, 
1979;  Bard  and  Bouchon,  1980a.b,  1985;  Harmsen  and  Harding,  1981;  Dravinski, 
1983;  Lee  and  Langston,  1983;  Sanchez-Sesma,  1983;  Nowack  and  Aki,  1984;  Bard 
and  Gariel  1986;  Geli  et  al.,  1987;  Koketsu,  1987;  McLaughlin  et  ah,  1987;  Moczo  et 
ah,  1987;  Mossessian  and  Dravinski.  19S7,  1992;  Benz  and  Smith,  1988;  Bravo  et  ah 
1988;  Campillo  et  ah,  1988;  Ivawase,  1988;  Sanchez-Sesma  et  ah,  1988;  Vidale  and 
Helmberger,  1988;  Ivawase  and  Aki.  1989;  Rial,  1989;  Seligman  et  ah,  1989;  Hill  et 
ah,  1990;  Horike  et  ah,  1990;  Novaro  et  ah,  1990;  Gaffet  and  Bouchon,  1991;  Gariel 
et  ah,  1991:  Koketsu  et  ah.  199U  Liu  et  ah,  1991:  Papageorgiou  and  Kim,  1991; 
Frankel  and  Vidale,  1992;  Graves  end  Clayton.  1992;  Kagawa  et  ah,  1992;  Ivawase 
and  Sato,  1992:  Ohori  et  ah.  1992;  Toshinawa  and  Ohmachi,  1992;  Uebayashi  et  ah, 
1992;  Frankel,  1993;  Jongmans  and  Campillo,  1993;  Mateos  et  ah,  1993;  Zahradnik 
et  ah,  1993).  Other  types  of  complex  geological  structures  investigated  include 
irregular  subsurface  layering  (Smith.  1975;  Kelly  et  ah,  1976;  Cerveny  et  ah,  1977; 
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Psencik,  1979;  Drake,  1980;  Nicc-letis,  1981;  Chapman  and  Drummond,  1982; 
Cerveny  and  Psencik,  1983,  1984;  Hill  and  Levander,  1984;  Mikhailenko,  1984,  1985; 
Virieux,  1984,  1986;  Levander  and  Hill,  1985;  Reshef  and  Ivosloff,  1985;  Dablain, 
1986;  Beydoun  and  Keho,  1987;  Campiilo  1987a;  George  et  al.,  1987;  Hu  et  al.,  1988; 
Paul  and  Campiilo,  1988;  Reshef  et  al.,  1988;  Bouchon  et  al.,  1989;  Chavez-Garcia 
and  Bard,  1989;  Daudt  et  al.,  1989;  Yamanaka  et  al.,  1989;  Carcione,  1993;  Gibson 
et  al.,  1993;  Schultz  and  Toksoz,  1993),  laterally-varying  crustal  structure  (Aki  and 
Larner,  1970;  Boore,  1970.  1972;  McMechan  and  Mooney,  1980;  Cerveny  et  al.,  1984; 
Muller.  1984;  Stephen,  1984:  Helmberger  et  al.,  1985a, b;  Vldale  et  al.,  1985; 
Campiilo,  1987b;  Emmerich,  1989,  1992;  Ho-Liu  and  Helmberger,  1989;  Maupin, 
1989;  Regan  and  Harkrider,  1989a, t;  Vaccari  et  al.,  1989;  Regan,  1991;  Campiilo  et 
al.,  1993  ;  Gregersen  and  Vaccari,  1993),  and  fault  zones  (Cormier  and  Spudich, 
1984;  Cormier  and  Beroza,  1987:  Ben-Zion  and  Aki,  1990;  Kang  and  McMechan, 
1993:  Moczo  and  Bard,  1993). 

In  the  present  paper  we  apply  the  boundary  integral  equation  method  to 
investigate  the  propagation  of  seismic  waves  over  distances  of  a  few  hundred 
kilometers  in  laterally  heterogeneous  crustal  structures.  We  use  a  formulation  very 
close  to  the  one  developed  by  Kawase  (1988)  and  Ivawase  and  Aki  (1989,  1990) 
where  the  Green's  functions  are  calculated  by  the  discrete  wavenumber  method  and 
where  the  singularities  inherent  to  the  boundary  integral  equation  approach  are 
removed  by  integrating  analytically  the  Green's  function  expressions  over  surface 
elements. 

Description  of  the  method 

The  simplest  medium  configuration  involves  two  homogeneous  media  separated 
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by  an  interface  and  is  depicted  in  Figure  la.  A  source  of  elastic  disturbance  is 
located  in  one  of  the  medium  (denoted  thereafter  by  index  1)  and  produces  at  a 
point  P  of  medium  1  a  direct  wavefield  V0(P )  .  For  simplicity  we  shall  only 
consider  the  two-dimensional  anti  plane  problem.  Using  Huygens  principle,  the 
wavefield  diffracted  by  the  interface  can  be  described  as  the  radiation  from 
secondary  sources  distributed  all  along  the  interface.  The  total  wavefield  at  P  may 
thus  be  written  in  the  form: 

V(P)-V0(P)  +  fo(Q)G(P,Q)dQ  (1) 

5 

where  o(Q)  is  a  source  density  function  which  represents  the  strength  of  the 
diffracting  source  at  the  interface  paint  Q  and  G(P,Q )  is  the  wavefield  radiated  at 
P  bv  a  unit  source  located  at  Q  snd  is  called  the  Green’s  function  of  medium  1. 
The  integration  is  taken  over  the  surface  of  separation  S. 

Similarly,  the  wavefield  diffracted  in  the  second  medium  may  be  written  in  the 
form: 

r(P')  =  /  o>[Q)G'(P',Q)dQ  (2) 

5 

where  G'  is  the  Green’s  function  of  medium  2. 

The  first  step  towards  the  obtention  of  a  numerical  solution  for  the  diffracted 
wavefield  requires  a  discretization  of  the  surface  integral.  This  is  achieved  by 
approximating  the  surface  5  by  N  surface  elements  AS,-  on  which  the  source  density 
fonctions  a  and  d  are  assumed  to  be  constant  (Figure  lb).  Equations  (l)  and  (2) 
thus  become: 

nn-  v„(p)+  i>,  /  o(p,Q)dQ  (3) 

i-l  AS, 
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nn  -  s  ^  /  a\p',Q)dQ 

i—l  AS, 

Then  choosing  P  and  P'  to  lie  0:1  the  interface  itself  and  denoting  by  Q:  this 
particular  point  (chosen  for  instance  at  the  middle  of  the  j th  surface  element, 
Figure  lc),  one  gets: 

V{Qj)  =  Vo(Qj)  +  I  G(Qj,Q)dQ  (4) 

i-l  AS, 

(4) 

V(Qj)=  /  G\Qj,Q)dQ 

i-l  AS, 

The  continuity  of  the  displacement  wavefield  across  the  interface  requires  the 
equality  of  the  two  right  hand  sides  of  equations  (4)  and  provides  a  system  of  N 
equations  (as  j=l,N)  where  the  unknowns  are  the  (7,-  and  the  </,-  .  The  continuity 
of  the  stresses  across  5  provides  N  more  equations  and  thus  leads  to  a  system  of  2 N 
equations  to  2 N  unknowns. 

Before  inverting  the  system  one  need  to  evaluate  the  expressions: 

U;j  -  /  G{Qj,Q)dQ  (5) 

45, 

and 

G'ij  =  /  G\QrQ)dQ 

45, 

for  t=l.Ar  :  ;=l.Ar 
and 
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(6) 


Tij  -  / 

AS, 


T'ij  =  /  T'(Qj,Q)dQ 

AS, 


■where  r(Q;fQ)  is  the  surface  traction  produced  at  Qj  by  a  unit  line  force  acting  at 
Q  ,and  V  is  the  same  quantity  calculated  for  the  elastic  parameters  of  medium  2. 

Using  the  discrete  wavenumber  method  (Bouchon  ana  Aki,  1977b),  the  two- 
dimensional  antiplane  Green  s  function  expressed  in  the  frequency  domain  may  be 
written  in  the  form: 


GiQjiQ)  = 


1  A/  £ 

2 ipprL 


with: 


k  _2tt 
r'  L  ’ 


Ur.  =  irj  —kn~)l/'2,  Im  (~tn  )<0  . 


where  (xj,zj)  are  the  coordinates  of  Qj  and  ( xQ,zQ )  those  of  Q,  p  is  the  medium 
density.  8  is  the  shear  wave  velocity.  L  is  the  periodicity  length  associated  with  the 
method,  and  M  is  an  integer  large  enough  to  insure  convergence  of  the  series. 

The  traction  T  is  thus  given  bv: 


dij  dZj 

v here  nx j  and  n. j  denote  the  i  ani  t  components  of  the  normal  to  the  interface  at 
Qj.  If  we  approximate  each  surface  element  AS,-  by  a  plane  surface,  the  analytical 
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integration  of  G  and  T  (equations  (5)  and  (6))  in  their  discrete  wavenumber  form  is 


straightforward  as  long  as  i=Fj.  When  i=j  equation  (8)  breaks  down  for  Q—Q^, 
and  equation  (6)  must  be  replaced  by  (see  appendix): 


(9) 


The  inversion  of  the  linear  system  expressing  the  continuity  of  displacement 
and  stress  across  the  interface  then  leads  to  the  two  source  density  distributions  cr(- 
and  c?i  representing  the  diffracted  wavefield  in  the  two  media. 


The  use  of  the  discrete  wavenumber  method  of  calculating  Green’s  functions  in 
the  boundary  element  scheme  has  the  advantage  of  preventing  the  occurrence  of 
mathematical  or  numerical  singularities  often  associated  with  boundary  element  or 
boundary  integral  equation  techniques.  It  is  also  easy  to  implement  when  one  of  the 
medium  or  the  two  media  are  made  up  of  flat  layers.  In  such  a  case,  the  full-space 
Green’s  functions  are  simply  replaced  by  the  layered  medium  Green’s  functions 
through  the  use  of  a  Thomson-Haskell  or  Kennett  algorithm  (Ivennett,  1974;  Muller, 
1985)  and  the  surface  integrals  (5)  and  (6)  are  still  evaluated  analytically.  We  shall 
now  present  such  examples  of  calculation. 

Test  of  accuracy  of  the  method 

In  order  to  check  the  precisioa  of  the  method  we  consider  the  configuration 
depicted  in  Figure  2a.  The  medium  consists  of  a  homogeneous  crustal  layer,  30  km 
thick,  overlaying  a  mantle  half-space.  The  source  is  a  line  of  horizontal  shear 
dislocation  occurring  on  a  vertical  plane  and  located  at  10  km  depth.  The  receivers 
are  placed  along  a  linear  profile  which  extends  in  a  direction  perpendicular  to  the 
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line  of  dislocation.  The  time  dependence  of  the  dislocation  is  a  smooth  step- 
function  defined  bv: 

f  (t)=[l-batan(t/t0  )J/2  with  a  rise  time  t0  equal  to  0.5s  (io) 

Two  calculations  are  made:  For  the  first  one  we  consider  the  problem  as  one 
involving  a  source  embedded  in  a  fiat  layered  medium  and  use  the  discrete 
wavenumber  method  coupled  with  refectivitv  matrices.  For  the  second  calculation 
we  consider  that  we  have  two  independent  layered  media  separated  by  a  fictitious 
vertical  interface  located  at  200  km  from  the  source.  We  thus  treat  the  problem  as 
if  the  crustal-mantle  structure  on  both  sides  of  the  200  km  mark  were  different.  We 
divide  the  fictitious  vertical  bounds  ry  into  100  surface  elements.  We  calculate  the 
mathematical  expressions  of  the  Green's  functions  G  and  G'  for  the  crust-mantle 
structure  using  the  discrete  wavenumber  method.  We  integrate  analytically  over 
each  surface  element  (equations  (5)  and  (6))  the  resulting  expressions.  We  invert  the 
resulting  linear  system  of  equations  and  obtain  the  two  source  distributions  o  and  o'. 
We  finally  use  the  two  source  distributions  inferred  to  calculate  the  seismic 
displacement  produced  at  the  receivers.  In  carrying  out  this  procedure  we  assume 
that  the  fictitious  surface  of  separation  between  the  two  media  extends  from  the 
free  surface  down  to  a  finite  depth  (chosen  as  45  km)  below  which  little  seismic 
energy  is  present. 

The  comparison  between  the  two  sets  of  results  is  displayed  in  Figure  2b.  The 
calculation  is  made  over  a  time  wirdow  of  60s  and  for  frequencies  up  to  2Hz.  The 
periodicity  length  L  used  in  the  discrete  wavenumber  method  for  the  two 
calculations  is  850  km.  The  agreement  between  the  two  solutions  proves  the 
validity  of  the  approach.  The  choice  of  the  number  of  elements  to  represent  the 
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diffracting  surface  (100)  is  somewhat  arbitrary.  In  general,  the  number  of  elements 
depends  on  the  particular  frequency  considered:  At  low  frequencies  a  minimum 
number  of  elements  is  required,  while  at  high  frequencies  this  number  should  be 
chosen  such  that  at  least  three  surface  elements  are  sampled  per  seismic  wavelength. 
For  this  reason  the  element  size  may  vary  from  one  layer  to  the  next  according  to 
the  layer  shear  wave  seismic  wavelength. 

Effect  of  a  vertical  fault 

We  now  investigate  how  the  propagation  of  Lg  waves  is  affected  by  the 
presence  of  a  fault.  Below  periods  of  a  few  seconds  and  up  to  very  high  frequencies 
the  Lg  wave  train  is  the  most  prominent  seismic  phase  produced  by  crustal 
earthquakes  at  regional  distances  (100km  to  1000km  from  the  source).  These  waves 
are  made  up  of  shear  waves  multipl'-reflected  in  the  crust  and  incident  on  the  Moho 
at  angles  more  grazing  than  the  critical  angle  defined  as  the  incident  angle  beyond 
which  all  the  downgoing  shear  energy  is  reflected  back  into  the  crust.  The 
prominence  of  the  Lg  waves  arises  from  the  waveguide  nature  of  their  propagation 
and  it  is  interesting  to  investigate  how  irregularities  of  the  w-aveguide  affect  their 
amplitude  and  characteristics. 

The  crustal  model  considered  i:;  showm  in  Figure  3a.  A  vertical  fault  with  2km 
offset  extends  from  the  surface  to  the  Moho.  The  fault  is  represented  by  boundary 
elements  and  the  Green's  functions  are  calculated  for  the  two  flat-layer  structures 
present  on  both  sides  of  the  fault.  The  criterion  used  to  determine  the  number  of 
elements  is  the  one  described  in  the  previous  section.  The  boundary  extension  is 
limited  to  45km  depth.  The  earthquake  is  modelled  as  a  line  of  horizontal  shear 
dislocation  occurring  on  a  vertical  olane  at  a  depth  of  10km  and  is  located  200km 
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from  the  fault.  The  time  dependence  of  the  dislocation  is  governed  by  equation  (10) 
and  the  frequency  range  extends  from  0  to  3Hz.  The  receivers  are  placed  along  a 
linear  array  perpendicular  to  the  line  of  dislocation  and  to  the  fault.  The 
corresponding  seismograms  are  disolaved  in  Figure  3b.  The  presence  of  a  back 
scattered  wavefield  originating  from  the  fault  is  clearly  seen.  Its  amplitude  is  about 
l/10th  the  amplitude  of  the  primary  field.  At  the  crossing  of  the  fault  a  change  in 
the  character  of  the  seismograms  occurs  with  the  near-diseappearence  of  the  surface 
wave  trains.  Also  a  shift  in  the  arrival  time  of  the  refracted  mantle  shear  wave  (Sn) 
takes  place  when  crossing  the  fault. 

A  comparison  between  the  seismogram  obtained  300km  away  from  the  source 
with  those  that  would  be  obtained  in  the  absence  of  the  fault  is  shown  in  Figure  4. 
The  first  comparison  (CRUSTl)  corresponds  to  the  case  of  a  flat-layer  crust  similar 
to  the  one  present  on  the  left  hand  side  of  the  fault.  The  second  one  (CRUST2) 
corresponds  to  the  crustal  model  on  the  right  hand  side  of  the  fault.  These  results 
show  that  the  presence  of  the  fault  changes  entirely  the  Lg  waveform  but  does  not 
affect  significantly  its  energy.  These  observations  can  be  accounted  for  by  the 
nature  of  the  Lg  phase.  Because  its  waveform  is  the  result  of  the  interference 
pattern  between  numerous  shear  wave  arrivals,  relatively  small  changes  in  crustal 
structure  produce  phase  differences  between  the  interfering  waves  and  thus  affect 
the  waveform.  On  the  other  hand  the  amplitude  of  each  individual  shear  wave 
constituting  the  Lg  wave  train  is  rot  significantly  affected  by  a  small  variation  in 
crustal  structure  so  that  the  total  energy  of  the  Lg  wave  is  almost  unchanged. 

Effect  of  a  change  in  the  Moho  depth 

We  now  investigate  the  effect  c  n  the  seismograms  of  a  lateral  change  in  crustal 
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thickness.  The  structure  considered  is  depicted  in  Figure  5a.  It  involves  a  change 
in  the  Moho  depth  which  rises  from  35km  to  25km  over  an  horizontal  distance  of 
10km.  The  earthquake  source  and  the  receivers  have  the  same  characteristics  and 
locations  as  in  the  previous  example.  The  position  of  the  diffracting  boundary  used 
in  the  calculation  is  represented  by  the  dashed  line.  It  extends  from  the  free  surface 
down  to  47km  depth  and  includes  the  dipping  Moho.  The  verticallity  of  the  two 
boundary  segments  extending  in  the  crust  and  in  the  mantle  is  chosen  to  minimize 
the  number  of  boundary  elements  required. 

The  results  are  presented  in  Figure  5b.  They  show  that  the  back  scattered 
wave  field  is  very  strong  which  suggests  that  in  real  observations  part  of  the  Lg 
wave  coda  is  due  to  back  scattering  enduced  by  variations  in  crustal  thickness  and 
that  the  level  of  coda  present  miglr;  be  related  to  the  roughness  of  the  Moho.  This 
simulation  also  shows  that  the  poshion  where  the  change  in  crustal  thickness  takes 
place  may  be  inferred  by  a  well-placed  seismic  array. 

The  comparison  of  the  seismogram  obtained  at  300km  with  those  that  would  be 
obtained  if  the  crust  were  laterally  homogeneous  is  depicted  in  Figure  6.  The  two 
flat-layer  models  correspond  to  the  structures  on  both  sides  of  the  region  where  the 
variation  of  crustal  depth  occurs.  Like  for  the  previous  configuration,  the  high 
sensitivity  of  the  Lg  waveform  to  the  crustal  structure  is  apparent,  while  the 
amplitude  stays  stable.  On  the  contrary  the  surface  waves,  which  sample  only  the 
upper  part  of  the  crust,  are  in  phase  between  the  three  models  and  display  similar 
waveforms.  The  arrival  time  of  the  refracted  mantle  shear  wave  (Sn)  lies  between 
the  arrival  times  predicted  by  the  two  fiat-layer  models. 

In  order  to  gain  more  physical  insight  into  the  diffraction  phenomenon  induced 
by  the  change  in  crustal  thickness,  we  present  in  Figure  7  a  space-time  display  of 
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the  wavefield.  The  amplitude  of  the  waves  is  calculated  at  a  grid  of  points 
extending  from  the  surface  down  to  a  depth  of  45km  and  between  epicentral 
distances  of  150km  and  250km.  Fifteen  snapshots  taken  at  3.5s  interval  and 
starting  39s  after  the  occurrence  of  the  earthquake  are  presented.  The  first  ones 
show  the  propagation  of  the  shear  wave  fronts  which  form  the  Lg  wave.  The  back 
propagation  of  the  scattered  travefi  ild  becomes  apparent  after  about  63s.  The  last 
five  snapshots  show  the  presence  of  a  nearly-cylindrical  wavefront  which  seems  to 
originate  from  the  free  surface  area  located  above  the  dipping  Moho.  We  interpret 
this  feature  as  shear  waves  scattered  upward  by  the  Moho  slope  and  subsequently 
reflected  downward  by  the  free  surface.  Because  of  the  dominance  of  the  primary 
wavefield  at  earlier  times,  only  the  downward  propagation  is  visible. 

Conclusion 

We  have  investigated  the  propagation  of  Lg  waves  in  a  laterally-varying  crust 
using  a  numerical  simulation  technique  formulated  in  terms  of  boundary  integral 
equations  where  the  Green’s  functions  are  evaluated  by  discrete  wavenumber 
summation.  The  method  is  well  suited  for  studying  the  propagation  of  seismic 
waves  in  a  layered  crustal  structure  which  varies  laterally  over  part  of  the 
propagation  path. 

We  have  considered  the  casejwhere  a  fault  with  vertical  offset  is  present  and 
the  case  of  a  change  in  crustal  thickness.  Our  results  confirm  the  simulation 
experiments  of  Maupin  (1989),  Regan  and  Harkrider  (1989a),  Campillo  et  al.  (1993), 
and  Gregersen  and  Vaccari  (1993)  which  all  clearly  show  that  Lg  wave  propagation 
in  crustal  structures  with  strong  lateral  variations  is  surprisingly  robust  and  that 
large  scale  geometric  features  of  the  crust  are  not  sufficient  to  account  for  the 
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strong  attenuation  of  Lg  waves  observed  in  many  regions. 

The  present  results  also  suggest  a  possible  relation  between  the  level  of  Lg  wave 
coda  and  the  degree  of  roughness  of  the  Moho.  They  further  indicate  the 
importance  of  the  back  scattered  wave  field  and  suggest  its  possible  use  to  map 
crustal  heterogeneities. 
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Appendix 


The  Somigliana  representation  theorem  (see  e.g.  Aki  and  Richards,  1980)  gives 
the  expression  of  the  displacement  or  stress  fields  as  a  function  of  their  value  on  an 
interface  S  delimiting  a  volume  V  .  In  the  SH  case  for  the  interior  problem,  the 
displacement  is  written  as: 

cuy(f)  =  /[<?w(x,f)rw(x)  -  ryw-(x?0«y(x)]d5y(x)  (la) 

s 

where  Gyy  and  Tyyj  are  the  displacement  and  stress  Green’s  functions  evaluated  at 
x  for  a  source  at  £  .  The  volumic  body  forces  are  assumed  to  be  null  and  c  is  a 
constant  which  takes  the  value  0,  1  or  0.5  for  x  respectively  outside,  inside  V ,  and 
on  S,  assuming  S  has  smooth  boundaries.  Different  equivalent  representations  may 
be  derived  from  (la)  to  obtain  a  form  equivalent  to  Huygens  principle  which  involves 
a  single  surface  distribution  (see  e.g.  Coutant,  1989).  In  the  SH  case,  the  Kirchhoff- 
Helmholtz  representation  for  instance  uses  a  single  layer  force  distribution  (or 
unidirectional  forces  as  opposed  to  a  double  layer  dipole  with  moment  distribution) 
whose  amplitude  is  given  by  a  stress  discontinuity  [t).  In  this  case  the  displacement 
is  written: 

«,(x)  =  /G„(x^)Mf/(9<lSy(9  <2a) 

S 

By  construction  of  this  representation,  the  displacement  field  is  continuous  across  5 
and  its  expression  is  valid  for  x  located  on  5  or  inside  or  outside  V.  The  stress  field 
representation  however  is  discontinuous  across  S,  Tout—rin  =  [rj,  and  its  value  on  S 
must  be  evaluated  with  a  limiting  process  (see  e.g.  Sanchez-Sesma  and  Campillo, 
1991).  The  stress  field  representation  is: 
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ryj  (x )  ?• 


.  M»(») 


+  /rntf(x^)[V(Qrfs*(0 

5 


where  f  is  1  when  x  is  on  S,  0  otherwise.  The  discretization  of  integral 
representations  (2a)  and  (3a)  is  achieved  by  approximating  the  surface  S  by  N 
planar  segments  AS,-  with  normal  n-,  and  by  assuming  *he  force  density  per  surface 
unit  <^=[i]yjnj  to  be  constant  on  each  segment.  The  representation  then  becomes: 


«,(*)  =  I  Gss(x,QdS(Z) 

«-l  AS, 


ry,(x)  - 


Mw(x)]  .  a 


+  /  Ty3(x,Z)dS(Z) 

i- 1  AS, 


In  order  to  solve  the  boundary  conditions  problem,  we  need  to  compute  the 
stress  and  displacement  on  the  boundary  S  at  the  N  segment  collocation  points  Xj. 
The  following  auto-influence  terms  have  to  be  evaluated: 


«»(*l)  "  <*i  /  Gyy(X|  ,0^(0 
A  S, 


rv(*i)  = 


<X5) 


+  ai  I  rw(xi,£)&?(0 


The  two  integrals  contain  respectively  a  weakly  and  a  strongly  singular  kernel 
and  can  be  integrated  in  the  space  domain  (e.g.  Sanchez-Sesma  and  Campillo,  1991) 
or  using  the  discrete  wavenumber  (DW)  decomposition  of  the  Green’s  function. 
Performing  the  integral  in  the  DW  domain  requires  however  a  minor  correction  to 
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the  expression  of  the  stress  field  as  a  function  of  the  horizontal  wavenumber.  The 
SH  displacement  Green’s  function  for  an  infinite  space  may  be  expressed  in  its 
wavenumber  representation  form  as: 


Gyyix)  —  .  / 

Anipp  _o, 


oo  -nf|2-z,|  -ik(x-X') 


- dk 


The  stress  tensor  Tyyy  horizontal  wavenumber  representation  can  be  directly 
obtained  from  (8a)  by  derivation  whh  one  exception,  when  the  source  and  "receiver" 
are  located  at  the  same  vertical  position,  ie  for  z~zg  .  In  this  case  the  component 
Tyyz  is  undefined  and  the  stress  can  only  be  computed  for  a  vertically  oriented 
(^=1,712=0)  segment.  Since  the  integral  is  independent  of  the  segment  orientation, 
we  evaluate  (7a)  by  performing  the  calculation  for  a  vertical  orientation,  that  is: 


f  Tyy(x\,£)dS(£)  -  /  Tm(x„ QdS(Q 

AS,  AS, 


1  °° 
4~ ip£P  AS.-oo 


■dkdz  —  0  (9a) 


This  result  is  similar  to  the  one  obtained  by  integration  in  the  space  domain,  and  we 
finally  get: 

(  \ 

’VW-T"  UOa) 
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Figure  legends 


Figure  1 

Illustration  of  the  method 
Figure  2 

Comparison  of  the  surface  displacement  traces  obtained  using  boundary  elements 
with  the  flat  layer  solution.  The  configuration  used  is  displayed  in  (a).  The  star 
indicates  the  source  location.  Receivers  are  located  at  epicentral  distances  between 
100  and  300km,  at  10km  interval.  The  two  superposed  sets  of  seismograms  are 
presented  in  (b).  A  reduced  time  equal  to  the  epicentral  distance  divided  by  the 
mantle  shear  wave  velocity  has  been  applied  to  the  traces. 

Figure  3 

Effect  of  the  presence  of  a  vertical  fault.  The  crustal  model  and  source-receivers 
configuration  are  depicted  in  (a).  The  surface  displacement  traces  are  presented  in 
(b).  A  reduced  time  equal  to  the  epicentral  distance  divided  by  the  mantle  shear 
wave  velocity  has  been  applied  to  tie  seismograms. 

Figure  4 

Comparison  of  the  seismograms  obtained  at  an  epicentral  distance  of  300km  in 
presence  of  a  fault  (model  of  Figure  3a)  with  those  obtained  for  the  two  fiat  layered 
structures  situated  on  the  left  hanc  side  (crust  1)  and  on  the  right  hand  side  (crust 
2)  of  the  fault. 

Figure  5 

Effect  of  a  change  in  crustal  thickness.  The  crustal  model  and  source-receivers 


46 


configuration  are  depicted  in  (a).  The  surface  displacement  seismograms  are 
displayed  in  (b)  with  a  reduced  time  based  on  the  mantle  shear  wave  velocity. 

Figure  6 

Comparison  of  the  seismograms  obtained  at  an  epicentral  distance  of  300km  for  the 
laterally-varying  crustal  model  of  Figure  5a  (Moho)  with  those  calculated  for  the 
crustal  model  with  constant  35km  thickness  (crust  3)  and  the  one  with  25km 
thickness  (crust  4). 

Figure  7 

Snapshots  showing  the  propagation  of  the  seismic  wave  field  for  the  configuration 
depicted  in  Figure  5a.  Each  frame  represents  a  cross-section  of  the  medium  between 
depths  of  0  and  45km  and  between  epicentral  distances  of  150  and  250km.  The  first 
snapshot  is  taken  39s  after  the  source  emission,  and  the  subsequent  ones  are 
displayed  at  3.5s  intervals. 
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FIGURE  1 
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FIGURE  2 
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II  -  Modeling  of  seismic  waves  generated  by  nuclear  tests  in  the  NTS  and 
Taourirt  Tan  Afella,  Hoggar  : 

II  - 1  Teleseismic  wave  form  modeling  including  geometrical  effects  of  superficial 
geological  structures  near  the  seismic  sources  at  the  NTS: 

We  analyse  observed  seismograms  of  21  events  recorded  at  teleseismic  distances  in 
France  from  nuclear  explosions  detonated  at  Nevada  Test  Site  (NTS).  Variations  of  the 
displacement  waveform,  duration,  and  amplitude  are  studied  in  terms  of  influence  of 
the  explosion's  medium  of  burial  and  location  in  the  laterally  heterogeneous  Yucca  Flat 
basin.  The  analysis  is  made  using  numerical  modelings  of  data  which  simulate  the 
modifications  of  the  original  source  signal  by  die  geometry  of  the  geological 
surrounding  structure.  Spalling  and  non-linear  effects  are  not  included  in  computations. 
The  numerical  simulations  of  the  variations  are  processed  using  a  mixed  symbolic  and 
numerical  algorithm  developed  to  simulate  the  ground  motions  that  may  be  recorded  in 
any  kind  of  two  dimensional  heterogeneous  non-vertical  structures.  This  algorithm  is 
linked  (/)  to  the  discrete  wavenumber  -  boundary  integral  equation  method  and  to  the 
reciprocity  theorem  to  simulate  the  displacement  field  radiated  to  farfield  distances  by 
the  source  site  geological  structure,  and  (//')  to  ray  propagation  to  propagate  the 
displacement  field  across  the  mantle  from  the  source  region  to  the  receiver.  Variations 
of  the  computed  displacement  amplitudes  are  as  large  as  twice  from  one  detonation 
point  to  another.  The  shapes  of  the  observed  seismograms  are  modelled  with  a  good 
reliability  for  most  of  the  21  events  in  the  frequency  band  studied,  i.e.  0.2H z  -  2.5Hz.  A 
set  of  2 1  relative  yield  estimates  are  derived,  which  include  the  source  site  response 
and  thus  the  amplitude  variations  induced  by  the  heterogeneous  structure  of  the  source 
region. 


II  -  2  Modeling  of  french  nuclear  tests  in  Taourirt  Tan  Afella  Massif,  Hoggar, 
Sahara : 

The  French  nuclear  test  site  for  underground  explosions  was  located  in  the  Hoggar 
Massif  at  the  beginning  of  the  1960's,  4  km  west  from  In  Eker,  and  150  km  north  from 
Tamanrasset  in  the  South  of  Algeria.  The  test  site  is  located  in  a  granitic  massif,  the 
Taourirt  tan  Afella  (5°2'  E;  24°3'  N),  which  is  intruded  in  two  gneiss  series.  The  massif 
is  intrusive  in  the  west  side  of  a  1  km  thick  mylonitic  corridor.  This  corridor 
approximately  oriented  north-south  separates  the  In  eker  gneiss  serie  (west  side)  and 
the  Tefedest  gneiss  serie  (east  side).  Both  series  present  a  north-south  schistosity. 
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Teleseismic  waveform  modeling  including  geometrical  effects 
of  superficial  geological  structures  near  to  seismic  sources 

Stephane  Gaffet 

Institut  de  Geodynamique  -  250.  rue  Albert  Einstein  -  Sophia- Antipolis  1  -  06560  Valbonne  -  France 
Abstract 

We  analyse  observed  seismograms  of  21  events  recorded  at  teleseismic  distances  in 
France  from  nuclear  explosions  detonated  at  Nevada  Test  Site  (NTS).  Variations  of 
the  displacement  waveform,  duration,  and  amplitude  are  studied  in  terms  of  influence 
of  the  explosion’s  medium  of  burial  and  location  in  the  laterally  heterogeneous  Yucca 
Flat  basin.  The  analysis  is  made  using  numerical  modelings  of  data  which  simulate  the 
modifications  of  the  original  source  signal  by  the  geometry  of  the  geological  surrounding 
structure.  Spalling  and  non-linear  effects  are  not  included  in  computations.  The  numerical 
simulations  of  the  variations  are  processed  using  a  mixed  symbolic  and  numerical 
algorithm  developed  to  simulate  the  ground  motions  that  may  be  recorded  in  any  kind 
of  two  dimensional  heterogeneous  non-vertical  structures.  This  algorithm  is  linked  (i) 
to  the  discrete  wavenumber  -  boundary  integral  equation  method  and  to  the  reciprocity 
theorem  to  simulate  the  displacement  field  radiated  to  farfield  distances  by  the  source  site 
geological  structure,  and  (ii)  to  ray  propagation  to  propagate  the  displacement  field  across 
the  mantle  from  the  source  region  to  the  receiver.  Variations  of  the  computed  displacement 
amplitudes  are  as  large  as  twice  from  one  detonation  point  to  another.  The  shapes  of  the 
observed  seismograms  are  modelled  with  a  good  reliability  for  most  of  the  21  events  in  the 
frequency  band  studied,  i.e.  0.2Hz  -  2.5Hz.  A  set  of  21  relative  yield  estimates  are  derived, 
which  include  the  source  site  response  and  thus  the  amplitude  variations  induced  by  the 
heterogeneous  structure  of  the  source  region. 

INTRODUCTION 

It  has  long  been  recognized  that  the  analysis  of  seismic  waveforms  induced  bv  under- 
ground  nuclear  tests  could  provide  valuable  information  for  the  understanding  of  source 
surrounding  site  effects.  Conversely,  these  source  site  effects  must  be  taken  into  account 
to  obtain  reliable  source  information  from  seismograms. 

Neither  the  influence  of  physical  properties  of  the  source  medium  (e.g.  density, 
porosity,  water  saturation)  that  must  be  taken  into  account  in  source  function  modeling  of 
underground  nuclear  tests  (Werth  and  Herbst,  1963;  Boardman  et  al.  ,  1964:  Hasegawa, 
1971,  1972;  Mueller  and  Murphy,  1971,  1972;  Power,  1974),  nor  non-linear  reflections 
or  spalling  effects  axe  studied  in  this  paper,  i.e.  we  focus  our  analysis  to  displacement 
waveform  variations  observed  at  teleseismic  distances  that  can  be  ascribed  to  scattering 
by  structure  heterogeneities  in  the  vicinity  of  the  source  region  (Cleary  et  al.  .  1975; 
Bouchon,  1976:  McLaughlin  et  al.  ,  1987;  Ferguson.  1988;  McLaughlin  and  Jih,  1988; 
Stead  and  Helmberger,  1988). 

Figure  1  displays  the  seismograms  recorded  by  the  French  Laboratoire  de  Detection 
et  de  Geophysique  (LDG)  network  for  detonations  in  the  central  region  of  Yucca  Flat, 
Nevada.  These  seismograms  were  constructed  by  delay  and  sum  of  5  vertical  recordings 
from  the  central  part  of  France.  The  epicentral  distance  is  79°-  81°  and  the  azimuth  from 
North  is  39°  -  41°  .  The  epicentral  distance  corresponds  to  an  angle  of  incidence  from  the 
source  region  of  16.95°  (Pho  and  Behe,  1972). 

Waveforms  differ  significantly  from  one  explosion  to  another.  For  instance,  Breton 
(top  in  Figure  1)  and  Texarkana  (bottom)  display  short  duration  signals  (2  seconds) 
while  Jornada  or  Dalhart  show  much  longer  durations  (more  than  6  seconds).  Since  the 
propagation  paths  from  the  source  region  to  the  receivers  are  identical,  and  since  the 
receiver  response  remains  constant  for  all  these  explosions,  the  perturbation  and  the 
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duration  lengthening  of  the  observed  waveforms  can  be  related  to  reflected  paths  and  to 
scattering  from  geological  structures  located  near  to  the  sources.  We  then  limit  our  study 
to  the  farfield  influence  of  the  medium  structure  surrounding  the  source  and  therefore 
ignore  the  receiver  structures  in  this  paper. 

In  a  first  part,  we  describe  a  new  mixed  symbolic  and  numerical  procedure  developped 
to  compute  teleseismic  synthetic  waveforms  for  sources  located  in  heterogeneous  media. 
The  numerical  procedure  is  an  extension  of  the  discrete  wavenumber  boundary  integral 
equation  method  (Bouchon  and  Aki,  1977).  For  this  extension,  any  kind  of  non-vertical 
subsurface  or  topographic  structures  can  be  introduced  in  the  model.  In  a  second  part, 
we  study  variations  of  the  amplitudes  and  the  waveform  perturbations  induced  by  the 
heterogeneous  structure  of  Yucca  Flat  as  seen  at  the  LDG  network.  We  show  that  our 
synthetic  seismograms  actually  reproduce  the  observed  variations  in  amplitude  as  well 
as  the  anomalous  durations  of  the  observed  seismograms  and  that  the  source  medium 
description  used  might  be  accurate  enough  to  describe  the  global  waveforms  recorded. 

NEAR  SOURCE  SITE  EFFECTS  EXPECTED  AT  YUCCA  FLAT 

The  map  of  Paleozoic  basement  depths  at  Yucca  Flat  (Figure  2  from  Ferguson,  1988)  shows 
a  North-South  elongated  basin  with  large  depth  variations  in  the  East- West  azimuth.  The 
locations  of  the  21  studied  explosions  are  depicted  by  circles.  A  geological  cross  section, 
AA’,  normal  to  the  main  axis  of  the  basin  and  roughly  turned  toward  the  France,  is  shown 
in  Figure  3.  At  this  cross  section,  the  width  of  the  valley  is  approximately  12.5km  and  the 
maximum  depth  to  Paleozoic  basement  is  1  km.  The  water  table,  wt  ,  is  from  Doty  and 
Thordarson  (1983).  It  separates  the  dry,  DT  ,  and  the  wet,  WT  ,  tuff  levels.  Above  these 
volcanic  rocks  are  superficial  deposits,  SD  ,  whose  base  depth,  sd  ,  is  interpolated  from 
drill  hole  data  given  by  Fernald  et  a l  (1968).  Following  Bouchon  (1976),  a  superficial 
layer  made  up  of  fan  alluvium,  FA  ,  is  introduced.  Finally,  the  Tertiary-Paleozoic  contact 
depth,  tp  ,  is  derived  from  the  works  of  Taylor  (1983),  Ferguson  et  al.  (1988)  and  Ferguson 
(198S).  The  cross  section  shows  the  complexity  of  the  geological  units  in  the  central  region 
of  Yucca  Flat.  The  aim  of  this  paper  is  to  investigate  whether  these  near  source  structure 
heterogeneities  could  account  for  the  observed  waveforms  recorded  in  France. 

A  basic  effect  of  near  source  heterogeneities  is  a  non-isotropic  radiation  pattern  of  the 
energy.  As  a  result  the  source  site  effect  influences  on  mb ,  can  be  cancelled  by  averaging 
records  from  stations  well  distributed  (in  azimuths  and  distances)  around  the  source.  But 
by  constrast,  if  only  records  from  narrow  aperture  network  are  available,  such  as  from  the 
LDG  network,  it  is  of  critical  importance  to  take  the  near  source  site  effects  into  account 
(McLaughlin  et  al.  ,  1987,  McLaughlin  and  Jih,  1988). 

MODELING:  METHOD  OF  COMPUTATION 

The  farfield  P- wave  displacement,  DIS]ar ,  is  assumed  to  result  from  the  convolution  of 
the  source  function,  SOU ,  with  the  source  site  response  SIT,0u,  with  the  propagation  ath 
response,  PRO.  with  the  receiver  site  response,  S7rree,  and  with  the  recording  svst«m 
response,  INS  .  Thus 

DISjar  =  SOU  *  PRO  *  INS  *  SITrtc  *  SIT,0v  (1) 

Different  formulations  of  the  source  function  SOU  have  been  proposed  in  terms  of  reduced 
displacement  potential,  (see  for  example  the  Haskell  model  (1967)  and  two  modified  forms 
by  von  Seggern  and  Blandford  (1972)  and  by  Helmberger  and  Hadley  (1981)).  The  source 
formulation  used  in  this  paper  is  the  one  given  by  Mueller  and  Murphy  (1971),  which 
takes  into  account  both  the  jdeld  and  the  detonation  depth.  This  model  involves  scaling 
factors  that  depend  on  the  source  medium  and  which  have  been  empirically  determined 
for  different  source  environments  (Viecelli,  1973;  Aki  et  al.  ,  1974;  Murphy,  1977;  Denny 
and  Johnson,  1991). 
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The  first  few  seconds  of  the  teleseismic  observations  at  distances  of  interest  (i.e.  SO0) 
essentially  consist  of  paths  that  propagated  across  the  mantle.  Assuming  ray  propagation 
in  the  mantle,  the  attenuation  along  these  paths  cam  be  described  by  the  Futterman  (1962) 
anelastic  attenuation  formulation,  which  in  the  frequency  domain,  is  given  by 

PRO  =  exp(-irfim  {1  -  erp-^))exp(-2tV/r  (Q  -  ^Logej-)Jexp^2ift’LogtfJ  (2) 

where  /  is  frequency,  t‘  =  Jpath  ^p-ds  represents  the  ray  path-integrated  effect  of  the  mantle 
quadity  factor  Q  ,  /„  is  the  frequency  above  which  dispersion  occurs  in  the  mantle  (f„  >  0), 
e  is  the  Euler  constant.  In  the  0.2Hz  -  2.5Hz  frequency  bamd  the  t‘a  value  can  vary  between 
0.1  and  0.8  (Frasier  and  Filson,  1972;  Cormier,  1982;  Der  and  Lees,  1985).  In  (2)  the  first 
exponential  scales  the  amplitude  according  to  the  factor  t“,  the  second  term  introduces 
a  phase  delay,  and  the  last  term  represents  the  frequency  dispersion  function  with  a 
phase  advance  for  /  >  lHz  and  a  phase  delay  for  /  <  lHz.  This  description  of  the  ray 
attenuation  in  the  mantle  implies  that  the  geometrical  spreading  remains  constant  for  all 
the  explosions  studied,  which  is  the  case  in  the  source-receiver  configuration  studied. 

The  response  of  the  LDG  network  short  period  recording  system,  INS ,  is  dominated  by 
the  seismometer  response  function  with  a  characteristic  frequency  of  lHz  and  a  damping 
factor  of  l/\/2- 

The  receiver  response,  SITre e,  remains  constant  for  all  the  explosions  studied  and  is 
assumed  to  be  flat  over  all  the  frequenc2r  band  for  the  next  computations.  The  term  SIT}OU , 
computed  using  the  discrete  wavenumber  -  boundary  integral  equation  method,  contains 
the  response  of  the  site  geological  structure  surrounding  the  source. 

Numerical  method  for  source  site  effects  modeling 

Different  numerical  methods  have  been  investigated  to  study  the  influences  of  the  near 
source  heterogeneities  on  the  teleseismic  P-wave  seismograms.  Hasegawa  (1971  and  1972) 
developed  a  method  and  correlated  the  recorded  seismogram  shape  to  the  geological 
complexity  of  the  source  environments.  Bouchon  (1976)  proposed  an  alternative  method 
which  is  related  to  the  Thomson-Haskell  method  (Thomson,  1950;  Haskell,  1953)  and  to 
the  reciprocity  theorem.  He  demonstrated  the  great  influence  of  the  source  depth  and  of 
the  burying  medium  on  the  mj  estimates.  The  main  limitation  of  these  methods  is  that 
the  medium  must  be  made  up  of  flat  layers.  Such  a  simplification  does  not  apply,  even 
approximately,  to  the  Yucca  Flat  basin  (Figure  2). 

Another  method  proposed  by  Aki  and  Larner  (1970)  has  been  used  by  Ferguson 
(1988)  for  Yucca  Flats  explosion  modelings,  and  other  empirical  or  data  analysis  methods 
(Taylor,  19S3:  Lay  and  Welc.  1987;  Lay  1987a  and  1987b)  emphasize  the  effect  of  near 
source  heterogeneous  structures  on  the  farfield  displacement  records.  These  methods  only 
use  a  few  frequencies  to  estimate  the  displacement  spectrum,  and  hence  to  estimate  the 
yield. 

The  algorithm  developed  below  removes  the  restrictions  of  the  previous  methods 
which  assume  flat-layered  media  or  are  limited  to  a  few  amplitude  samples  computed 
in  the  frequency  domain.  This  algorithm,  in  association  with  the  discrete  wavenumber  - 
boundary  integral  equation  method  and  with  the  reciprocity  theorem,  is  able  to  compute 
the  complete  waveform  in  the  time  domain  for  events  located  in  any  kind  of  topography 
and  subsurface  structure. 

Continuity  of  the  stress-displacement  field  across  the  interfaces 

We  assume  that  the  structure  of  the  medium  studied  is  entirely  described  by  a  set  of 
interfaces  as  shown  in  Figure  3.  An  interface  is  defined  as  the  boundary  that  separates 
two  given  media  (one  of  the  two  media  may  be  the  vacuum,  i.e.  the  boundary  is  a  free 
surface).  Each  medium  is  thus  generallv  bounded  by  more  than  one  interface.  We  use  the 
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matrix  formulation  to  describe  the  problem  and  do  not  make  explicit  the  content  of  these 
matrices.  The  reader  may  refer  to  Gaffet  and  Bouchon  (1989),  Bouchon  et  al.  (1989)  and 
Gaffet  and  Bouchon  (1991)  for  the  complete  mathematical  formulation. 

In  the  discrete  wavenumber  -  boundary  integral  equation  method,  each  side  of  an 
interface  i  ,  is  associated  with  a  distribution  of  regularly  spaced  horizontal  and  vertical 
forces  (see  Figure  4).  The  notation  Q{ ,  is  used  to  represent  the  vectorial  forces  located  on 
each  side  of  the  interface  i  in  contact  with  the  media  I.  The  potentials  radiated  by  these 
forces  are  used  to  compute  the  field  diffracted  by  all  of  the  interfaces.  The  field  diffracted 
from  an  interface  j  by  its  forces  located  in  the  medium  J  and  received  by  the  side  of  am 
interface  i  in  contact  with  the  medium  I  is  given  by  the  product  A^Qj .  At  this  first  step 
of  the  description,  the  interface  i  can  receive  the  diffracted  field  in  medium  I  from  the 
interface  j  if  and  only  if  {iff )  I  -  J-  The  dimension  of  the  matrix  Aj,-  equals  4.V7,  x  2 M, 
in  the  P  -  SV  problem  and  2A/,-  x  Mj  in  the  SH  problem,  M,  and  Mj  being  the  number  of 
points  used  to  define  the  interfaces  i  and  j. 

We  introduce,  now,  a  term  written  Sf  that  describes  a  seismic  source  located  within 
the  medium  7,  that  acts  on  the  interface  i.  Using  the  notations  previously  defined,  we  can 
write  the  equation  which  represents  the  interactions  between  the  interface  i  and  all  the 
interfaces  that  describe  the  model. 


5/  +  A fiQ-  +  £  Aj:Qj~r  =  S*  +  Aii  Qi  +  £  Aji Qj  /  >  *7  *  separates  2  media :  (4) 


S(  +  A (iQi  +  Y2  A jiQjSl  =  0.  if  I  ln  (4)  is  vacuum 


(5) 


Initiation  of  the  algorithm:  first  step  of  elimination  procedure 

The  equations  (4)  and  (5)  are  repeated  for  each  interface  i  and  each  product  of  type  A 
is  computed  iff  I  —  J.  (4)  and  (5)  are  rewritten  below  to  describe  the  system  of  matrix 
that  initiate  the  first  step  of  the  algorithm,  i.t.  level  /  =  0.  The  following  series  of  matrices 
and  vectors  is  obtained  for  each  interface  i 


°a,  is  a  4 Mi  x  4 Mi  matrix  (i.e.  P-SV  case)  or  a  2 Mi  x  2 Mi  matrix  (z.e.  SH  case)  that  is 
constructed  by  concatenation,  as  indicated  by  the  vertical  bar,  of  the  two  rectangular 
matrices  A(t  and  -Af, .  Then,  for  each  interface  j  acting  on  the  interface  i  .  we  have  to 
compute 


°b//  =  “a"1  [4  |  ~Aji 


(U 


The  terms  Ajt-  exist  iff  it  has  been  computed  in  (4)  and  (5),  i.e..  iff  at  least  one,  non 
vacuum  ,  of  the  two  media  J  or  J  .  in  contact  with  interface  j  is  also  in  contact  with  i. 
Thus  for  each  interface  i,  we  obtain 


Qi  =  °ft  +  £  °b//  Qj  (8) 

jJjSi 

In  (8),  the  vectors  Qi  and  Q;  correspond  to  the  forces  applied  on  each  sides  of  the 
interfaces  i  and  j. 
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Symbolic  elimination 

We  solve  the  lineax  system  of  matrices  described  in  (8)  using  symbolic  Gauss  Jordan 
elimination.  The  unknowns  are  the  vectors  Q,  and  the  elements  of  the  system  are  the 
rectangular  matrices  bjf.  The  system  may  be  sparse  depending  on  interaction  between 
interfaces  i  and  j ,  t.e.  depending  on  the  topology  of  the  structure.  Using  (8),  at  step  l  of 
elimination,  we  obtain 

Qi  =  '-1A  +  Q,  4-  Y,  ,_1  b'/  <?i.  f°r  *’  >  1  (10) 

3>iJ& 

The  main  difficulty  is  to  avoid  multiplication  between  null  terms  at  step  /,  i.e.  the  terms 
corresponding  to  non  interactive  interfaces.  An  automatic  symbolic  procedure  is  adopted 
which  considers  that  the  product  ,_1bjjJ  Qi  must  be  taken  into  account  iff  there  exists  a 
fictive  path,  called  Ci,  that  may  cross  the  interfaces  1  to  /  to  connect  the  interfaces  /  and 
i,  allowing  the  diffracted  field  emitted  by  the  interface  l  in  the  medium  L  or  L'  to  reach 
the  interface  i  in  the  medium  I  or  i' .  The  same  condition  is  applied  to  take  into  account 
the  product  l~lbjf  Qj.  A  path,  called  C2,  must  cross  the  interfaces  1  to  /  to  allow  to  the 
field  diffracted  by  the  interface  j  in  the  medium  J  or  J  to  be  received  by  the  interface  i 
in  the  medium  I  or 

This  mixing  of  simultaneously  numerical  calculation  and  symbolic  determination  whether 
a  computation  has  to  be  made  is  the  crux  of  the  proposed  method. 

In  a  similar  way  we  also  express  Qi  as 

Qi  =  '"‘A  +  '-'bff  Qi  +  £  '^b/f  Qj  (11) 

The  product  l~lhffQi  must  now  be  taken  into  account  iff  a  path,  called  Cj,  exists  and 
crosses  the  interfaces  1  to  l  which  permits  at  step  l  to  the  field  diffracted  by  the  interface 
i  in  the  medium  I  or  /'  to  be  received  by  the  interface  l  in  the  medium  L  or  L‘ .  Finally,  to 
take  into  account  the  product  ,-1b j,LQj,  a  path,  called  C3,  must  exist  across  the  interfaces 
1  to  l  which  allows  the  field  diffracted  by  the  interface  j  in  the  medium  J  or  f  to  be 
received  by  the  interface  l  in  medium  L  or  if . 

Then  if  we  replace  Q\  in  (8)  by  its  expression  given  in  (9)  we  get 

Qi  =  (  l~l0i  +  '_1b^  '-1#  )  +  ‘-'bi;1  '-'bff  Qi 

+  E  <  +  '‘lb"  '"‘b“  )  (12) 

and  thus 

«,  =  [/-  ’-‘bf'  '-b',1  ]"‘( '-v,  +  '-b£'  '-■« ) 

+  £  [  /  -  ’-■b*'  '-‘b?  ]’l(  ‘-‘b//  +  >-‘b V  '-‘b jf- )  Q,  (is) 

j>lj& 

We  then  obtain 

Qi  =  l0i  +  'b//  Q;  (14) 

i>ij& 
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The  different  products  in  (13)  that  define  the  matrices  'b//  and  the  vectors  :Q,  in  (14)  are 
computed  iff  the  corresponding  paths  Ci,  Cj,  C3,  or  C3  can  be  found,  i.e.  ,  we  determine 
symbolically  whether  or  not  the  products  of  matrices  in  (13)  have  to  be  numerically 
calculated.  We  have  described  the  iteration  /  -  1  to  /  of  the  matrix  linear  system  reduction 
process  that  has  now  to  be  repeated  for  /  =  1  onto  the  Nth  interface  to  complete  the 
resolution  of  (8). 

Backpropagation  of  the  forces 

At  the  end  of  this  elimination  process,  i.e.  (10),  ....  (14),  we  obtain  an  upper  triangular 
system  of  matrices.  Its  resolution  is  straightforward  that  we  call  backpropagation  of  the 
forces.  The  last  iteration  in  (14)  gives  directly  the  force  vector  QNtk..  This  vector  is  then 
introduced  in  the  Nth  -  1  expression  to  compute  and  so  on  to  back  compute  all 

the  forces  that  must  be  applied  on  both  sides  of  each  interface  in  order  to  insure  the 
continuity  of  the  stress-displacement  field  across  the  interfaces  considered. 

Note  that  this  method  described  here  may  actually  be  applied  to  solve  any  kind  of 
linear  band  system  with  the  minimum  need  of  computer  memory  (each  matrix  may  be 
saved  on  mass  storage  until  it  is  needed),  and  with  the  maximum  of  numerical  stability. 

NUMERICAL  MODELING  OF  YUCCA  FLAT  DATA 

The  map  given  in  Figure  2  displays  the  distribution  of  the  21  studied  detonations.  Figure 
3  shows  the  corresponding  cross  section  A  A'  used  for  this  studv.  The  elastic  parameters 
of  each  medium  summarized  in  Table  1  are  those  from  Ferguson  (1988).  o,  3,  p  and  v 
denote  the  P  and  SV  velocities,  the  density  and  Poisson’s  ratio. 

The  focal  parameters  of  the  explosions  are  listed  in  Table  2  (from  USGS  bulletin). 
The  locations  of  the  21  explosions  are  projected  onto  the  vertical  cross  section  AA’  in 
Figure  5.  We  assume  that  the  21  locations  have  accurate  depths  and  therefore  determine 
the  medium  of  burial.  Under  this  assumption  4  explosions  are  detonated  above  thw  water 
table  in  tuff  (e.g.  Draughts  ,  Tajo  ,  Strake  .  and  Texarkana  ).  Most  of  them  are  below  the 
water  table.  Breton  is  assumed  to  be  buried  in  the  sediment  deposit  level  ( SD  ,  Figure 
1)  and  appears  to  be  very  near  the  closing  area  (C.4  )  of  all  the  media  near  the  horst 
Paleozoic  structure. 

Influence  of  the  source  location 

We  compute  the  seismograms  that  should  be  recorded  in  France  for  a  line  of  explosions 
buried  at  a  depth  of  650m  in  the  wet  tuff  layer,  and  for  yields  of  lkt,  lOkt,  and  lOOkt.  The 
Mueller  and  Murphy  (1971)  source  function  is  used  with  the  parameters  corresponding 
to  the  wet  tuff-rhvolite  medium  (Murphy,  1977).  The  t~  value  is  taken  to  be  0.7  (Cormier, 
1982).  The  maximum  amplitudes  and  the  location  of  the  corresponding  explosions  are 
presented  in  Figure  6.  We  obtain  a  strong  variation  by  a  factor  of  2  of  the  maximum  0  to 
peak  amplitude  with  the  location  of  the  detonations  (Figure  6  bottom),  which  appears 
to  be  insensitive  to  the  energy  of  the  explosions  {i.e.  for  lkt,  lOkt.  and  lOOkt  plotted 
respectively  as  black  losanges,  plus  signs,  and  crosses). 

The  amplification  is  maximum  for  detonations  located  in  the  western  zone  ( WZ  Figure 
6).  In  this  area,  the  amplitude  variation  is  very  sensitive  to  the  location  of  the  detonations: 
amplification  can  vary  by  as  much  as  40%  for  a  horizontal  shift  of  only  a  few  hundred 
meters  of  the  detonation  point.  The  minimum  amplification  is  observed  for  explosions 
located  in  the  eastern  saturated  ( EZ  )  tuff  zone.  This  minimum  is  not  only  related  to 
the  shape  of  the  Paleozoic  basement  but  rather  appears  to  integrate  the  geometry  of  all 
the  surrounding  heterogeneous  structures.  The  maximum  attenuation  occurred  for  shots 
detonated  in  the  western  part  of  this  zone.  We  obtain  at  this  point  a  relative  attenuation 
larger  than  50%  if  compared  to  the  absolute  basin  maximum  amplification,  and  40%  when 
compared  to  the  eastern  zone  maximum  amplitude. 
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These  preliminary  results  confirm  the  conclusion  of  previous  investigation  regarding 
the  necessity  to  take  into  account  the  heterogeneous  structures  that  may  exist  in  the 
vicinity  of  seismic  sources  to  accurately  analyse  seismogram  amplitudes  obtained  from  a 
narrow  aperture  network. 

Waveform  fitting  -  Yield  inversion 

Synthetic  seismograms  have  been  computed  using  the  Mueller  and  Murphy  (1971)  source 
function  and  compared  to  our  records.  The  constraint  used  to  construct  the  synthetic 
waveforms  is  simply  to  best  fit  the  relative  amplitudes  of  the  observed  seismograms  by 
only  adjusting  the  yield.  These  best  fitting  seismograms  are  shown  in  Figure  7.  Rummy 
(m»  =  5.7),  which  has  the  maximum  amplitude  observed  was  chosen  as  the  reference 
event.  We  therefore  arbitrarily  fix  its  yield  to  be  loOkt,  This  calibration  then  allows  us  to 
estimate  the  yields  of  the  other  explosions.  The  two  last  columns  of  Table  3  summarize 
the  yields  obtained  using  this  calibration  and  the  observed  relative  amplitudes.  These  two 
columns  illustrate  the  yield  variations  that  may  occur  between  explosions  such  as  Hermosa 
and  Lowball  or  Rummy  and  Jornada  (detonated  in  wet  tuff),  or  Tajo  and  Texarkana 
(detonated  in  dry  tuff),  having  the  same  seismogram  amplitudes. 

The  quality  of  the  fits  is  quantified  using  the  correlation  coefficient  computed  over  the 
first  5  seconds  of  signal  (see  last  column  in  Figure  7).  Most  part  of  the  seismograms  have 
synthetic  waveforms  which  closely  follow  the  observed  ones.  The  longer  duration  of  the 
ground  motion  observed  for  events  in  the  central  region  (i.e.  around  Tahoka  )  is  correctly 
reproduced,  as  well  as  the  relative  amplitudes  of  the  two  first  troughs  corresponding 
to  the  P  and  pP  wave  arrivals.  The  fourth  peak  which  induces  duration  lengthening, 
arrives  approximately  4  seconds  after  the  first  arrival.  Such  a  delay  is  observed  on  the  real 
seismograms. 

There  are  different  ways  to  explain  the  slight  discrepancy  between  synthetic  pP-P 
and  observed  pP-P,  e.g.  a  too  shallow  depth  of  burial,  too  large  wave  velocities  or  lateral 
velocity  variations  as  reported  by  Johnson  and  McEvillv  (1990).  We  however  prefer  other 
explanations  which  concern  the  representation  of  the  source  radiation  for  an  explosion 
that  takes  place  near  an  interface  ( i.e .  at  distances  smaller  than  the  elastic  radius). 

(i)  For  numerical  modeling  assuming  linear  elasticity,  the  source  used  to  represent 
the  nuclear  explosion  must  not  be  considered  as  a  point  source  because  the  elastic  zone 
si2e  becomes  important  in  comparison  with  the  heterogeneity  dimensions.  Likewise,  the 
spherical  elastic  zone  shape  in  a  homogeneous  space  should  be  modified  in  presence  of 
different  media  having  different  compaction  properties  (Murphy,  1980).  Therefore,  we 
must  take  into  account  the  elastic  zone  shape  and  the  size  of  an  explosion  in  heterogeneous 
medium  to  explain  variations  in  the  frequency  content  and  in  the  amount  of  seismic 
energy  radiated  outside  of  the  basin,  (ii)  Another  reason  for  waveform  discrepancy  is  that 
secondary  sources  such  as  spall  and  non  linear  reflection  at  free  surface  are  not  included 
in  this  work  which  adresses  only  the  geometrical  influence  of  the  structure  surrounding 
the  source. 

Thus  if  the  shot  location  in  a  heterogeneous  medium  plays  a  great  role  in  shaping  the 
seismograms,  it  appears  necessary  to  take  care  in  the  explosion  source  model  to  have  the 
best  possible  source  as  input  into  the  diffracting  heterogeneous  media  system. 

CONCLUSIONS 

We  have  implemented  in  a  new  numerical  algorithm,  the  discrete  wavenumber  -  boundary 
integral  equation  method  to  simulate  source  site  effects  which  could  account  for  the  farfield 
radial  ground  motions.  The  computation  uses  the  reciprocity  theorem  to  simulate  21 
underground  explosions  from  the  central  region  of  Yucca  Flat  valley,  Nevada,  and  recorded 
in  France.  This  study  confirms  that  teleseismic  seismogram  waveforms  and  amplitudes  are 
very  sensitive  to  the  shot  location  for  explosions  located  in  a  heterogeneous  sedimentary- 
filled  basin.  We  have  obtained  a  good  correlation  between  both  shape  and  duration  of 
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synthetic  seismograms  and  the  observed  ones.  This  indicates  that  the  location  of  the 
explosions  is  accurately  taken  into  account  in  our  calculations  and  that  the  heterogeneities 
in  the  source  region  must  be  taken  into  account  for  waveform  analysis  and  for  yield 
estimates.  Source  site  effects  such  as  these  are  of  great  importance  in  calibrating  a 
receiver  response  in  relation  to  a  specified  site,  as  well  as  in  calibrating  the  site  of  a  single 
explosion.  We  show  that  if  the  near  source  heterogeneities  can  be  reliably  accounted  for, 
the  uncertainty  on  the  yield  estimate  of  a  nuclear  test  by  a  narrow  aperture  network  would 
be  significantly  reduced. 
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FIGURE  CAPTIONS 

Figure  1  -  Seismograms  recorded  in  France  by  the  French  L.  D.  G.  network.  The  time 
window  is  20  seconds  long.  Dates,  names,  and  relative  amplitudes  of  all  events  are  given 
above  each  seismogram. 

Figure  2  -  Map  of  Paleozoic  basement  at  Yucca  Flat  test  site,  Nevada.  The  circles  show 
epicenters  of  the  explosions  used  in  the  study.  The  two  first  letters  of  the  explosion  names 
are  also  given  (Table  1).  The  cross  section  A  A’  is  depicted  by  the  thick  bold  line.  This 
map  is  adapted  from  the  map  given  by  Ferguson  (19S8). 

Figure  3  -  Cross  section  AA’  (see  Figure  2  for  location).  Lower  case  letters  refer  to 
interfaces  (i.e.  sd  as  basement  of  sedimentary  deposits,  wt  as  water  table,  and  tp  as 
tertiary-Paleozoic  contact  depth).  Upper  case  letters  refer  to  geological  units  (i.e.  FA  : 
fan  alluvium,  SD  :  sedimentary  deposits,  DT :  dry  tuff,  WT :  wet  tuff,  and  PZ  :  paleozoic 
basement).  The  display  vertical  exaggeration  is  5:1  . 

Figure  4  -  Description  of  the  interfaces  and  media  notations  used  in  text.  We  show  here 
three  parts  of  three  interfaces  t  ,  j  ,  and  l .  The  couples  of  orthogonal  arrows  represent 
the  couples  of  horizontal  and  vertical  forces  applied  at  the  sampling  points  on  each  side 
of  all  interfaces.  The  name  of  the  media  separated  by  the  three  interfaces  are  I,  J,  f , 
L,  and  L' .  The  notations  used  in  text  to  represent  the  forces  applied  (i.e.  Qt,  Q;,  and  Qt) 
are  illustrated  in  this  Figure. 

Figure  5  -  Locations  of  the  21  explosions  studied  in  relation  to  the  cross  section  given  in 
Figure  3.  The  depths  and  locations  are  those  given  by  the  USGS  bulletin. 

Figure  6  -  Estimated  maximum  displacement  amplitudes  that  would  be  received  by  the 
French  LDG  network  for  a  set  of  18  fictitious  explosions  at  Yucca  Flat  test  site,  Nevada. 

Figure  7  -  Comparison  of  synthetic  seismograms  (thick  bold  line)  with  the  observed 
data  (thin  bold  line)  for  the  21  explosions  studied.  The  synthetics  are  10  seconds  long 
and  include  the  LDG  short  period  seismometer  response.  Dates  of  each  event  are  noted 
to  the  left.  The  bold  numbers  are  relative  amplitudes  recorded  in  France.  The  explosion 
names  are  followed  in  the  last  column  by  the  correlation  coefficient  between  synthetic  and 
observed  seismograms. 
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TABLE  1  -  Elastic  parameters  of  media 


Medium 

a(m/s)  3(m/s] 

I  p(kg/m3) 

V 

FA 

700 

330 

1500 

0.35 

SD 

1340 

640 

1770 

0.35 

DT 

2140 

1140 

1790 

0.30 

WT 

3000 

1600 

1790 

0.30 

PZ 

4570 

2640 

2500 

0.25 

TABLE  2  -  Event  locations  and  magnitudes 


No 

Latitude 

Longitude 

1 

37.0867 

-116.0711 

2 

37.0858 

-116.0686 

3 

37.0831 

-116.0661 

4 

37.0797 

-116.0514 

5 

37.0722 

-116.0500 

6 

37.0533 

-116.0453 

7 

37.0914 

-116.0511 

8 

37.0658 

-116.0464 

9 

37.0608 

-116.0453 

10 

37.0889 

-116.0492 

11 

37.0728 

-116.0458 

12 

37.0786 

-116.0439 

13 

37.0944 

-116.0447 

14 

37.0947 

-116.0361 

15 

37.0692 

-116.0303 

16 

37.0947 

-116.0322 

17 

37.0947 

-116.0278 

18 

37.0739 

-116.0197 

19 

37.0983 

-116.0156 

20 

37.0867 

-116.0069 

21 

37.0767 

-116.0006 

depth  (m)  Date  Event  name 

483  13  Sep  1984  5.0  Breton 

500  30  Aug  1988  5.0  Bullfrog 

600  22  Mar  1986  5.1  Glencoe 

442  27  Sep  1978  5.0  Draughts 

701  09  Nov  1977  5.7  Sandreef 

600  05  Dec  1985  5.7  Kinibito 

640  28  Jan  1982  5.9  Jornada 

639  01  Mar  1984  5.9  Tortugas 

600  13  Aug  1987  5-9  Tahoka 

600  13  Oct  1988  5.9  Dalhart 

533  14  Apr  1983  5.7  Turquoise 

564  12  Jul  1978  5.5  Lowball 

564  -25  May  1977  5.3  Crewline 

637  03  Jun  1975  5.7  Mizzen 

640  04  Feb  1976  5.8  Keelson 

640  02  Apr  1985  5.7  Hermosa 

594  27  Apr  1977  5.4  Bulkhead 

640  27  Sep  1978  5.7  Rummy 

500  05  Jun  1986  5.3  Tajo 

518  04  Aug  1977  5.0  Strake 

500  10  Feb  1989  5.2  Texarkana 
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TABLE  3  -  Yields  estimated 


No  depth  (m) 

m4  Event  name  buried 

in  V'0(kt)  A  /  .4, 

1 

483 

5.0  Breton 

SD 

57 

0.13 

2 

500 

5.0  Bullfrog 

WT 

28 

0.17 

3 

600 

5.1  Glencoe 

WT 

35 

0.23 

4 

442 

5.0  Draughts 

DT 

124 

0.24 

5 

701 

5.7  Sandreef 

WT 

139 

0.74 

6 

600 

5.7  •  Kinibito 

WT 

79 

0.43 

7 

640 

5.9  Jornada 

WT 

178 

0.98 

8 

639 

5.9  Tortugas 

WT 

135 

0.76 

9 

600 

5.9  Tahoka 

WT 

143 

0.77 

10 

600 

5.9  Dalhart 

WT 

170 

0.92 

11 

533 

5.7  Turquoise 

WT 

108 

0.59 

12 

564 

5.5  Lowball 

WT 

99 

0.54 

13 

564 

5.3  Crewline 

WT 

56 

0.31 

14 

637 

5.7  Mizzen 

WT 

93 

0.63 

15 

640 

5.8  Keelson 

WT 

87 

0.59 

16 

640 

5.7  Hermosa 

WT 

76 

0.54 

17 

594 

5.4  Bulkhead 

WT 

58 

0.41 

18 

640 

5.7  Rummy 

WT 

150 

1 

19 

500 

5.3  Tajo 

DT 

72 

0.20 

20 

518 

5.0  Stroke 

DT 

48 

0.15 

21 

500 

5.2  Texarkana 

DT 

52 

0.20 
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Yucca  Flat  cross-section  and  events  selected 
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Modeling  of  French  nuclear  tests 
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SUMMARY 

The  influence  of  the  topography  on  the  shape  and  amplitude  of  seismograms  recorded  at 
short  distances  is  investigated  for  a  set  of  nuclear  explosions  detonated  in  the  Touaxirt 
tan  Afella  mountain,  Hoggar,  Sahara,  mountain  where  was  located  the  French  nuclear 
test  site  during  1960’s.  The  wavefield  observed  in  one  side  in  comparison  to  the  other 
side  of  mountain  show  a  phase  generated  with  a  strong  amplitude  by  reflection  of  the 
direct  field  inside  of  the  mountain.  This  phase  is  back-diffracted  and  can  be  seen  in  North 
azimuth  and  not  in  South  azimut.  The  existence  of  such  a  phase  is  correlated  to  the  source 
location  inside  of  the  mounts'  i.  The  amplitudes  are  also  modified  by  the  topography  with 
a  variation  of  25%  between  observations  made  in  the  North  and  in  the  South  azimuths. 

Key  words:  nuclear  explosion,  numerical  modeling,  site  effect 

1  INTRODUCTION 

The  French  nuclear  test  site  used  at  the  beginning  of  1960’s  was  located  in  the  Hoggar 
massif  (Figure  1),  4  km  west  from  In  Eker,  and  150  km  North  from  Tamanrasset  in  South 
of  Algeria.  The  test  site  is  located  in  a  granitic  massif,  the  Taourirt  tan  Afella  (  5®2*E; 
24°3'N),  which  is  intruded  in  two  gneiss  series.  (Faure,  1972;  Boullier  and  Bertrand,  1981). 
The  massif  is  intrusive  in  the  west  side  of  a  lkm  thick  mylonitic  corridor  (Figure  2,  from 
Faure,  1972).  This  corridor,  approximately  oriented  North-South  separates  the  In  Eker 
gneiss  serie  (west  side)  and  the  Tefedest  gneiss  serie  (east  side).  Both  series  present  a 
North- South  schistosity.  The  massif  is  embedded  in  a  thin  level  of  sand  (Figure  3,  from 
Faure,  1972).  It  presents  an  ellipsoidal  shape,  8km  long  in  the  North- South  direction  and 
5.6km  large  in  the  east-west  direction  and  culminates  at  2000m  high  over  a  region  with 
1000m  mean  elevation. 

The  study  presented  here,  corresponds  to  the  first  step  realized  in  waveform  analysis 
of  the  French  explosions  detonated  in  1960’s  in  the  Taourirt  tan  Afella  massif,  Hoggar, 
Sahara.  This  study  concerns  the  influence  of  the  topography  onto  the  ground  motions  at 
local  distances  (  i.e.  from  lkm  to  30km)  in  the  aim  to  understand  the  waveforms  observed 
in  the  NS  and  EW  azimuths  at  BRI  (approximately  2km  southwest  from  the  explosions 
studied),  BRIII  (l5km  west),  INA  (  i.e.  In  Amguel,  30km  South),  and  by  the  French  LDG 
network  (2000km  in  the  North). 

2  DATA  COLLECTED 

13  underground  explosions  have  been  detonated  inside  of  the  massif  between  november 
1961  and  february  1966  (see  table  below).  Three  sets  of  the  ground  velocities  recorded,  are 
displayed  figures  4,  5,  and  6,  which  correspond  to  the  explosions  Rubis  ,  Opale  ,  and  Jade 
respectively.  Only  Rubis  has  been  recorded  at  short  distance  (l.73km)  from  the  source  at 
BRI .  All  explosions  have  been  recorded  at  distances  of  around  15km  at  BRIII  and  30km 
at  INA  . 

Place  for  Table  1 

The  shape  of  the  ground  velocity  differs  greatly  between  these  three  explosions.  The 
recordings  at  BRIII  and  INA  for  Rubis  show  pure  P  -wave  and  Rayleigh  wave.  The 
recordings  of  Opale  show  a  long  surface  wavetrain  onto  the  horizontal  component.  The 
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greatest  difference  is  found  for  recordings  of  explosion  Jade  .  It  presents  a  very  long 
and  strong  coda  associated  to  the  surface  wave  at  BRHI  and  IN  A  .  These  differences 
may  be  explained  by  different  causes,  i.e.  explosion  working  (spalling,  coupling),  water 
saturation  in  source  region  (and  its  evolution  between  each  experiment),  surrounding 
geological  structures  and  topography  which  induce  multipathing  and  amplification  effects. 

3  STUDY 

Our  interest  is  to  understand  the  ground  displacement  waveforms  at  local  distances,  lower 
than  30km,  in  both  NS  and  EW  azimuths,  and  at  teleseismic  distances  at  20s,  in  relation 
with  the  location  of  the  shot  inside  of  the  Taourirt  tan  Afella  massif,  and  with  the  azimuth 
of  recording.  In  the  aim  to  realize  this  objective,  the  Taourirt  tan  Afella  structure  is 
considered  as  to  be  superposition  of  (i)  the  topography  and  (ii)  the  crustal  underground 
structure  as  defined  by  Munier  (1982).  Assuming  linear  superposition  of  both  effects, 
this  decomposition  is  useful  to  separate  effects  induced  by  surficial  and  by  underground 
heterogeneities. 

The  results  presented  below,  are  relative  to  topographical  effects  of  Touarirt  tan 
Afella  for  the  NS  azimuth.  The}’'  discuss  about  the  ground  displacement  modifications 
induced  (i)  by  the  topography  and  (ii)  by  the  detonation  location  inside  of  the  mountain. 
The  numerical  method  used  for  computation  is  the  discrete  wavenumber  -  boundary 
integral  equation  method  (Bouchon  and  Aki,  1977;  Gaifet  and  Bouchon,  1989).  The  wave 
propagation  velocities  are  a  =  6km/s  and  0  =  o/\/3,  and  the  source  time  function  is  a 
Ricker  pulse. 

Figure  7  is  shown  as  reference,  the  ground  velocities  obtained  for  an  explosion  buried 
at  300m  depth  under  a  flat  surface.  The  only  waves  generated  are  the  direct  P  -wave 
followed  by  the  Rayleigh  wave  with  a  shorter  amplitude. 

Figure  8,  displays  the  ground  velocities  obtained  when  the  explosion  is  at  the  same 
depth  (  i.e.  300m),  straight  under  the  top  of  Taourirt  tan  Afella  mountain  (configuration 
is  given  figure  9).  The  source  pulse  is  centered  at  3Hz.  The  time  duration  is  13s  and  the 
observations  ranges  from  2lkm  in  South  (+  signs)  to  30km  in  North  (-  signs)  directions. 
Distances  and  relatives  amplitudes  are  written  for  each  H  and  Z  velocity  components. 
Besides  the  direct  field  (line  P,  Figure  8)  which  is  followed  by  a  Rayleigh  wave  (Rl) 
with  a  strong  amplitude  on  both  side  of  the  mountain,  two  other  branches  appear  ■which 
propagate  with  a  delay  of  1.7s  in  the  North  direction.  These  two  phases  are  (i)  reflexion 
inside  of  the  massif  of  the  direct  P  -wave  (RD)  and  (ii)  a  second  Rayleigh  wavetrain 
(R2)  generated  by  the  interaction  between  the  reflected  P  -wave  inside  of  the  irregular 
topography.  In  South  azimuth  a  phase  (S)  propagates  with  the  5  -wave  velocity.  It  is 
vertically  polarized  and  is  observed  at  the  front  of  the  Rayleigh  wave  branche  (Rl). 
This  phase  is  a  5  surface  wave  generated  in  the  flat  zone  of  the  massif.  The  maximum 
amplitude^  computed  in  North  direction  are  60%  the  one  obtained  in  the  South. 

Figure  10  shows  a  comparison  between  horizontal  (H)  and  vertical  (Z)  ground 
velocities  obtained  at  20km  North  and  South  from  the  source  in  the  configuration  shown 
figure  9.  Three  observations  can  be  made  using  this  explosion  location  configuration. 
(i)  The  amplitude  of  the  vertical  component  is  significantly  modified  by  the  topography- 
while  the  horizontal  one  is  not,  (ii)  the  shape  and  amplitude  of  the  direct  P  -wave  are 
not  affected  by  the  topography,  and  (Hi)  the  amplitude  and  shape  of  the  diffracted  field  ( 
i.e  Rayleigh  wave  and  secondary  phases)  of  both  horizontal  and  vertical  components  are 
different. 

We  now  study  the  influence  of  the  location  onto  the  ground  velocity  when  an  explosion 
is  detonated  in  South  (explosion  A,  figure  11)  and  in  North  (explosion  B,  figure  12)  inside 
of  the  mountain.  Explosions  A  and  B  are  supposed  to  be  buried  300m  depth  normal  to 
surface  topography.  The  distances  of  observation  is  approximately  15km  South  and  North 
from  the  source.  The  pulse  source  is  centered  at  4Hz,  and  the  seismograms  are  computed 
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for  a  time  window  of  7.8s  (Figure  13).  The  relative  amplitude  and  the  horizontal  offset  if 
written  in  the  right  side  of  each  component.  The  diffracted  field,  i.e.  mainly  composed 
by  the  phases  noted  (S),  (Rl),  (R2),  has  a  higher  frequency  content  for  explosion  A  than 
for  explosion  B.  This  is  observed  for  all  ground  velocity  components  on  both  sides  of  the 
massif.  This  observation  may  be  related  to  the  lower  elevation  of  the  mountain  in  the 
explosion  A  zone.  The  phase  (R'2)  is  well  observed  for  the  two  explosions  on  the  ArS  and 
BzN  components  while  the  surface  5  -wave  appears  on  the  A:N  and  BzS  components. 
This  behaviour  is  coherent  with  the  observation  previously  made,  figure  8.  We  notice  the 
phase  delay  of  the  (Rl)  phase  on  BhS  component. 

Thus,  the  different  phases  are  symetrically  observed  for  explosions  A  and  B.  But  as 
for  it,  the  frequency  contents  of  the  diffracted  fields  are  significantly  different,  between 
explosion  A  and  B,  as  well  as  the  relative  amplitude  of  the  different  components. 

4  CONCLUSION 

We  realize  here  the  first  step  of  a  global  seismological  study  of  Saharian  explosions.  The 
preliminary  results  show  the  great  influence  of  the  topography  onto  the  ground  velocity 
recorded  in  the  frequency  band  between  1  to  4Hz,  at  distances  from  1  to  30km  far  from 
explosions  detonated  inside  of  the  granitic  massif  of  the  Taourirt  tan  Afella,  Hoggar, 
Sahara.  '  '  ' 

The  general  result  is  that,  alone  the  topography  strongly  shapes  the  ground  velocities 
and  their  amplitudes.  This  result  encourages  us  to  enter  the  second  step  of  the  study 
concerning  the  underground  basaltic  and  mohorovicic  structures,  associated  to  the  lateral 
variations  of  the  medium,  such  as  the  NS  mylonitic  corridor,  in  the  aim  to  get  a  complete 
explanation  of  the  waveform  recorded. 

Yet,  the  main  results  are  the  following  ones.  The  back  scaterred  wave  shapes  very 
strongly  the  seismograms.  The  surface  wave  appears  to  be  generated  after  reflexion  of  the 
direct  field  onto  the  opposit  side  of  the  massif.  The  back  scaterred  field  is  made  up  of 
both  P  -surface  wave  and  Rayleigh  wave  which  simply  duplicate  the  direct  corresponding 
fields.  We  show  that  for  a  receiver  located  at  the  same  epicentral  distance,  the  amplitude 
of  the  vertical  ground  velocity  may  vary  with  a  factor  of  2,  from  one  side  to  the  other  side 
of  the  massif.  We  show  that  the  duration  of  the  surface  wavetrain  is  strongly  enlarged 
depending  on  the  side  (the  azimuth)  of  observation. 
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FIGURE  CAPTIONS 

Figure  1  :  Map  of  North  Africa.  The  gray  filled  circle  indicates  the  location  of  the  French 
nuclear  test  site  in  Algeria,  during  1960’s. 

Figure  2  :  This  map,  from  Faure  (1972),  illustrates  the  geological  formation  located 
around  the  Taourirt  tan  Afella  test  site. 

Figure  3  :  Topographical  section  in  North-South  azimuth  across  Taourirt  tan  Afella 
massif,  from  Faure  (1972). 

Figure  4  :  Seismograms  recorded  during  Rubis  experiment,  at  BRI ,  BRIII ,  and  IN  A  . 
The  time  scale  is  repeated  with  different  shape  for  each  trace. 

Figure  5  :  Seismograms  recorded  during  Opale  experiment,  at  BRIII  and  IN  A  . 

Figure  6  :  Seismograms  recorded  during  Jade  experiment,  at  BRIII  and  IN  A  . 

Figure  7  :  Horizontal  (H)  and  vertical  (Z)  components  of  the  ground  velocity  computed 
for  an  explosion  buried  300  m  depth  under  a  flat  surface  in  a  homogeneous  space.  The 
source  is  a  Ricker  pulse  centered  at  3.5  Hz.  The  time  window  is  7.8  s.  Receivers  are 
identified  by  their  epicentral  distances  (from  -18  km  to  +18  km). 

Figure  8  :  Horizontal  (H)  and  vertical  (Z)  seismograms  computed  in  the  configuration 
described  figure  9.  The  time  duration  is  13  s.  The  source  is  an  explosion  with  a  pulse  time 
function  centered  at  3  Hz.  The  epicentral  distances  of  the  receivers  are  from  +21  km  in 
the  South  to  -30  km  in  the  North.  The  relative  amplitude  are  written  above  the  epicentral 
distances.  The  notation  used  to  highlight  the  different  phases  are  discussed  in  the  text. 

Figure  9  :  Geometrical  configuration  used  to  compute  the  seismograms  presented  figure 
8.  The  P  and  S  wavelength  are  illustrated  above  the  topography. 

Figure  10  :  Comparison  of  ground  velocity  computed  at  an  epicentral  distance  of  20  km 
for  each  side  of  the  mountain.  See  text  for  explanation. 

Figure  11  :  Geometrical  configuration  used  for  the  explosion  A. 

Figure  12  :  Geometrical  configuration  used  for  the  explosion  B. 

Figure  13  :  Comparison  of  the  horizontal  and  vertical  components  of  the  ground  velocity 
for  the  explosion  A  (left  side)  and  for  explosion  B  (right  side).  For  explosion  A,  the  two 
receivers  are  14-9  km  far  in  the  South  and  in  the  North  from  explosion  location.  For 
explosion  B,  the  two  receivers  have  an  offset  of  1 5.5  km  in  the  North  and  in  the  South 
azimuths  from  explosion.  The  different  components  are  noted  as  follow  in  the  text,  AB 
hz  NS  .  The  first  letter  refers  to  the  explosion  (A  or  B),  the  second  letter  refers  to  the 
component  (H  or  Z),  and  the  third  letter  refers  to  the  azimuth  of  observation  (South  or 
North). 
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TABLE  1  -  French  explosions  characteristics 


Explosion 

Date 

Time  TU 

Lon  E 

Lat 

N  mb 

Agate 

11/7/1961 

llh  29mn  59.931s 

5°  3' 

7.6” 

24°  3' 

25.5" 

Beryl 

5/1/1962 

lOh  OOmn  0.458s 

5°  2 

30.8" 

24°  3' 

46.8" 

Emeraude 

3/18/1963 

lOh  02mn  0.351s 

5°  3' 

7.9" 

24°  2 

28.9" 

Amethyste 

3/30/1963 

9h  59mn  0.328s 

5°  3' 

25.2" 

24°  2 

36. 0" 

Rubis 

10/20/1963 

13h  OOmn  0.011s 

5°  2 

19. 0" 

24°  2 

7.8"  5.6  ( 

USGS) 

Opale 

2/14/1964 

llh  OOmn  0.347s 

5°  3' 

n 

8.6 

24®  3' 

13. l" 

Topaze 

6/15/1964 

13h  40mn  0.367s 

5°  2 

4.4" 

24®  3' 

59.8" 

Turquoise 

11/28/1964 

lOh  30mn  0.035s 

5°  2 

30. l" 

24®  2' 

30.7" 

Saphir 

2/27/1965 

llh  30mn  0.039s 

5®  1 

52.3" 

24®  3' 

31.4"  5.8  ( 

USGS 

Jade 

5/30/1965 

llh  OOmn  0.037s 

5°  3' 

3.l" 

24®  3' 

is.o" 

Corindon 

10/1/1965 

lOh  OOmn  0.043s 

5°  2 

2.6" 

24®  3' 

53.7 

Tourmaline 

12/1/1965 

lOh  30mn  0.088s- 

5°  i 

48.9" 

24®  2' 

37.4"  5.1  ( 

USGS ) 

Grenat 

2/16/1966 

llh  OOmn  0.035s 

5°  2' 

28.4" 

24®  2 

39.0" 

Figure  1 
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FIGURE  2 


FIGURE  5 


OP  ALE  (2/14/1964) 
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FIGURE  6 


JADE  (5/30/1965) 
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FIGURE  9 


a  velocity:  6.0  km/s  -  P  velocity:  3.5  km/s. 
Periodicity:  J00  km 
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Ground  velocity  for  a  pulse  source  centered  at  3  Hz  -  Time  window:  13  s 


FIGURE  10 


Ground  velocity  for  a  pulse  source  centered  at  3  Hz  -  Time  window:  13  s 
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FIGURE  11 


The  massif  is  embedded  in  a  thin  level  of  sand.  It  presents  an  ellipsoidal  shape,  8  km 
long  in  the  north-south  direction  and  5.6km  large  in  the  east-west  direction  and 
culminates  at  2000  m  high  over  a  region  with  1000  m  mean  elevation. 

The  study  we  have  carried  out  corresponds  to  the  first  step  in  wave  form  analysis  of  the 
French  explosions  detonated  in  1960's  in  the  Taouarirt  tan  Afella  Massif,  Hoggar, 
Sahara. 

The  study  is  composed  of  two  parts  : 

First  it  concerns  the  influence  of  topography  onto  the  ground  motions  at  local  distances 
(i.e.  from  1km  to  30km)  and  at  teleseismic  distances  (2000km).  Second,  it  concerns  the 
influence  of  the  underground  crustal  structures  at  the  same  distances. 

Both  steps  have  been  proceeded  in  the  aim  to  understand  the  waveforms  observed  in 
the  NS  and  EW  azimuths  at  station  BRI  (approximately  2km  Southwest  from  the 
explosions),  BRIII  (15  km  west),  INA  (i.e.  In  Amguel,  30  km  south),  and  by  the  French 
LDG  network  located  2000km  to  the  North. 

The  preliminary  results  show  the  large  influence  of  the  topography  onto  the  ground 
velocity  which  was  recorded  in  the  frequency  band  between  1  and  4  Hz,  at  distances 
from  1  to  30  km  far  from  explosions  detonated  inside  the  granitic  massif. 

The  general  result  is  that  the  topography  alone  strongly  shapes  the  ground  velocities 
and  their  amplitudes,  this  results  encourages  us  to  enter  the  second  step  of  the  study 
concerning  the  underground  basaltic  and  Mohorovicic  structures,  associated  to  the 
lateral  variations  of  the  medium,  such  as  the  NS  mylonitic  corridor,  in  the  aim  to  get  a 
comprehensive  explanation  of  the  recorded  waveform. 

The  main  results  are  the  following  : 

-  The  back  scattered  wave  shapes  very  strongly  the  seismograms.  The  surface  wave 
appears  to  be  generated  after  reflexion  of  the  direct  field  onto  the  opposite  side  of  the 
massif.  The  back  scattered  field  is  made  of  both  P-surface  wave  and  Rayleigh  wave 
which  simply  duplicate  the  direct  corresponding  fields. 

-  We  show  that  for  a  receiver  located  at  the  same  epicentral  distance,  the  amplitude  of 
the  vertical  ground  velocity  may  vary  with  a  factor  of  2,  from  one  side  to  the  other  side 
of  the  massif. 
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Ill  -  Data  processings  associated  with  a  mini-array  recently  built  in  the 
Center  of  France : 

III  - 1  Detection  and  phase  identification  capabilities  : 

During  the  80's,  substantial  efforts  have  been  carried  out  to  use  local  mini-arrays  for 
automatic  event  detection  (e.g.  :  Mikkelveit  et  al.,  1983).  Beside  these  studies,  some 
researches  have  also  been  undertaken  to  evaluate  their  capabilities  of  automatic 
azimuth  and  slowness  determination,  for  location  purpose. 

For  similar  objectives,  the  French  Laboratoire  de  Detection  et  de  Geophysique  (LDG) 
has  installed  in  1990  a  small  temporary  local  network  in  the  Center  of  France, 
provisory  composed  of  5  vertical  component  short  period  seismometers  with  an 
aperture  of  1.2  km  (In  the  next  future  10  stations  would  be  set  up).  Ninety  eight  (98) 
teleseismic  events  have  been  recorded  by  the  network  during  6  months  of  operation. 
We  present  here  the  main  result  concerning  the  automatic  determination  of  azimuth  and 
slowness  for  each  event  of  this  dataset.  Two  different  methods  of  data  processing  are 
tested  and  compared  for  that  purpose. 

1)  The  methods : 

•  The  first  method,  so-called  "frequency-wavenumber  (f-k)  method",  has  often  been 
used  by  seismologists  (e.g.  :  Capon,  1969,  Gupta  et  al.,  1990).  At  the  opposite  of  the 
original  use  of  this  method  which  computes  the  K-spectrum  using  filtered  signals  in  a 
narrow  band,  we  have  followed  the  algorithm  proposed  by  Nawab  et  al.,  (1985),  which 
uses  the  zero-delay  spectrum  to  obtain  a  k-spectrum  containing  information  integrated 
over  the  whole  frequency  range.  Then  the  azimuth  is  evaluated  by  searching  the 
maximum  value  of  the  radial  energy  of  the  K-spectrum. 

•  The  second  method  is  just  a  "correlation  method".  In  a  first  step,  the  cross-correlation 
function  of  the  stations  taken  two  by  two  leads  to  the  determination  of  arrival-time 
differences,  with  an  accuracy  of  one  sampling  interval  (i.e.  :  0.02  s  in  our  case).  In  a 
second  step,  the  cross-spectrum  phase  allows  to  compute  the  residual  arrival-time 
difference,  less  than  0.02  s.  This  residual  arrival-time  difference  is  determined  by  the 
slope  of  the  phase  as  a  frequency  function  in  the  characteristic  frequency  band  of  the 
signals. 
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Then,  assuming  a  plane  wave  as  a  propagation  model,  these  arrival-time  differences  are 
used  to  compute  both  azimuth  and  slowness  of  die  wave. 

2)  The  results  : 

We  have  processed  the  whole  dataset  of  98  teleseismic  events. 

Besides  12  events  (given  by  numbers),  the  standard  deviation  of  the  residuals  is  15 
degrees,  in  the  1000-8500  km  range.  For  larger  distances,  the  incident  wave  being  very 
close  to  the  vertical  axis  leads  to  a  poor  azimuthal  determination,  as  expected. 

A  similar  study  has  been  made  with  the  correlation  method,  for  the  determination  of 
azimuth  using,  first,  arrival-time  differences  obtained  from  the  cross-correlation 
functions,  and  secondly,  those  computed  by  both  cross-correlation  functions  and  cross¬ 
spectrum  phases.  It  clearly  shows  the  gain  obtained  by  the  use  of  the  cross-spectrum, 
specially  in  the  1000-8500  km  range.  This  is  again  demonstrated  by  testing  the 
consistency  of  the  arrival-time  differences  set  which  must  verify  the  triangular  Chasles 
relation  is  :  A tjj  -  Alik  +  A tkj  for  all  i;  j;  k.  The  RMS  value  of  the  residuals  of  the 
Chasles  relation  is  0.013  s  when  we  only  use  arrival-time  differences  determined  by 
cross-correlation,  and  0.002  s  in  the  second  case. 

Another  advantage  can  be  attributed  to  the  "  correlation  method".  At  the  opposite  of 
"(fk)  method"  which  assumes  a  plane  wave  and  model,  we  can  use  a  more  refined 
model  defined  by  a  plane  wave  and  a  set  of  time  delays  affecting  the  arrival-times. 
Using  the  global  dataset,  we  can  statistically  compute  each  station  anomaly  as  the  mean 
value  of  the  residuals.  This  leads  to  time  delays  ranging  from  -0.008  s  to  0.005  s 
producing  time  differences  greater  than  half  a  sampling  interval  for  some  couples  of 
stations. 

The  residual  azimuths  as  a  function  of  the  true  azimuths  obtained  from  USGS,  take  into 
account  the  station  anomalies,  clearly  shows  a  cosine  dependence  which  might  be 
explained  by  a  deeping  structure  of  the  crust  layers  below  die  network.  Final  results 
taking  into  account  this  cosine  dependence  within  the  range  of  1000-8500  km,  show 
that  the  standard  deviation  becomes  less  than  10  degrees  with  only  two  aberrant  points. 

We  have  used  a  set  of  98  telesismic  events  recorded  by  a  temporary  5  stations  mini¬ 
array  set  up  in  the  Center  of  France  to  test  two  methods  for  automatic  measurement  of 
P-wave  azimuth. 
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The  first  method,  or  "(f-k)  method",  gives  as  expected  consistent  results  in  most  of  the 
cases  (azimuth  determination  with  an  error  of  15  degrees  and  12  aberrant  evaluations). 

A  second  method,  the  "correlation  method",  derived  from  the  doublets  method 
(Poupinet  et  al.,  1982,  Plantet  et  al.;  1985),  uses  a  step-by-step  algorithm,  computing 
first,  on  set  time  differences,  then  uses  the  cross-spectrum  as  a  vernier  to  refme  these 
differences.  It  also  allows  to  introduce  station  anomalies  to  correct  the  propagation 
model  before  the  azimuth  computation.  Azimuth  is  then  determined  with  an  error  of 
less  than  10  degrees  in  the  1000-8500  km  range  with  only  two  aberrant  determinations. 
Consequently,  better  results  are  obtained  with  this  method. 

Further  studies  will  investigate  the  regional  domain  for  which  the  higher  frequency 
content  will  give  a  better  accuracy  in  the  arrival-time  differences. 

We  suggest  such  a  correlation  method  to  be  tested  with  seismic  data  recorded  in  other 
mini-arrays. 


Ill  -  2  Automatic  processing  of  seismic  events  recorded  on  a  mini-array  (signal 
analysis  combined  with  neural  networks)  : 

We  present  a  new  method  for  automatic  processing  of  seismic  events  recorded  on  a  5  - 
station  mini-array  located  in  Central  France. 

The  first  step  of  the  process  consists  in  the  computation  of  accurate  arrival  time 
differences  for  each  couple  of  stations  using  their  correlation  function.  These  arrival 
time  differencies  computed  for  different  time-windows  and  for  different  frequency 
bands  allow  us  to  get  three  time-frequency  plots  representing  the  consistancy  of  the 
time-delay  set,  the  velocity  and  the  azimuth  deduced  from  the  location  performed  with 
these  time-delays.  The  consistancy  of  the  time-delay  set  is  then  used  as  a  signal 
detector. 

For  teleseismic  events,  the  location  is  strickly  deduced  from  the  velocity  and  the 
azimuth  and  leads  to  an  accuracy  of  less  than  10  degrees  for  distances  lower  than  80 
degrees. 
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For  regional  events,  an  additional  step  is  needed  to  identity  the  different  phases.  This 
task  is  performed  by  a  neural  network  using  as  inputs  for  each  time-step  the  velocities 
computed  on  the  current  part  of  signal  filtered  in  different  frequency  bands  provided 
that  the  consistancy  of  the  time-delay  set  is  verified.  A  multi-layered  perceptron  then 
computes  the  possibility  of  appearance  for  each  phase  as  a  function  of  time.  The  more 
probable  azimuth  is  determined  and  used  as  a  filter.  Finally,  the  distance  is  computed 
by  choosing  in  propagation  tables  the  best  solution  according  to  the  possibilities. 
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Abstract.  Two  different  methods  devoted  to  automatic 
location  of  both  regional  and  teleseismic  events  recorded 
on  a  mini-arTay  are  compared  in  terms  of  location  accu¬ 
racy.  The  first  method  belongs  to  the  frequency- 
wavenumber  methods  family;  the  second  one  uses  signal 
processing  to  compute  accurate  time-delays  and  then  derive 
the  event  parameters.  In  the  latter  method  a  careful  study 
of  each  time-delay  set  is  performed  in  order  to  remove 
ambiguity  errors  using  Chasles  relationship.  It  leads  to  an 
accuracy  of  about  10  degrees  on  the  azimuth  determination 
for  distances  smaller  than  80  degrees,  and  3  degrees  for 
regional  distances.  Finally,  a  criterion  derived  from  these 
relations  is  proposed  for  phase  detection  and  identification. 

Introduction 


comparison  of  signals  recorded  by  different  stations  with  a 
high  accuracy. 

The  106  teleseismic  events  used  for  the  study  are  those 
which  have  been  localized  by  USGS.  The  epicentrai  dis¬ 
tance  range  extends  from  9  to  176  degrees,  and  twenty- 
eight  events  have  propagated  PKP  phases.  Their  magni¬ 
tudes  range  from  mb=4.1  (at  a  distance  of  10  degrees)  to 
mb=6.4. 

The  dataset  for  local  and  regional  events  is  composed 
of  30  earthquakes  recorded  and  localized  by  the  national 
permanent  LDG  network.  Their  epicentrai  distances  range 
from  42  to  1540  km  and  their  magnitude  from  ML=2.2  to 
Ml=5.3.  In  this  dataset,  for  16  events  the  Pg-phase  is  the 
first  arrival:  for  the  14  remaining  events,  Pn  is  the  first 
arrival. 


During  the  80*s,  many  attempts  were  made  to  adapt 
seismic  arrays  to  automatic  signal  detection  and  event  loca¬ 
tion.  Besides  these  studies,  some  research  has  been  con¬ 
ducted  in  order  to  evaluate  their  capabilities  of  automatic- 
azimuth  and  slowness  determination  [e.g.:  Mvkkeltveit  et 
al..  1990]. 

In  1989,  the  ’Laboratoire  de  Detection  et  de 
Geophysique’  (LDG)  temporarily  installed  a  5-station  small 
aperture  mini-array,  which  recorded  more  than  400  seismic 
events  during  its  4-month  operation  period.  This  dataset 
has  been  used  to  study  an  automatic  location  process 
adapted  for  both  teleseismic  and  regional  events.  Two 
methods  were  used  with  a  different  approach:  one  belongs 
to  the  spectral  methods  family,  the  other  analyzes  directly 
the  time  series. 

We  present  the  results  concerning  the  automatic  deter¬ 
mination  of  azimuth  and  slowness  for  two  sub-datasets 
(teleseismic  and  regional  events),  using  these  two  different 
methods,  and  compare  them. 

The  Data 

During  the  test,  the  mini -array  was  mainly  composed  of 
5  short-period  vertical  component  seismometers.  Figure  1 
shows  a  map  of  the  mini-array  centered  on  45.741 1°N 
2.0133°E.  It  has  an  extension  of  2.1  km  in  the  north-south 
direction  and  1.3  km  in  the  east- west  direction. 

The  recorded  signals  were  transmitted  to  station  GM5 
for  numerical  recording  after  a  12-bit  digitalization  at  a 
rate  of  50  samples  per  second.  This  structure  with  a  single 
clock  provides  an  homogeneous  digitalization,  and  allows 


The  Methods 

The  first  method  (denoted  by  'k-spectrum  Method-) 
belongs  to  the  ’frequency-wavenumber  methods’  family. 
The  classical  formuiation  of  these  methods  has  often  been 
used  by  seismologists  [e.g.:  Capon,  1969,  Ingate  et  al.. 
1985].  It  assumes  that  the  propagation  of  the  wave  through 
the  mini-arTay  can  be  modeled  by  a  plane-wave.  So  it 
investigates  the  horizontal  wavenumber  plane  (kxjcy)  for 
each  frequency,  and  determines  the  azimuth  and  the  slow¬ 
ness  of  the  wave  from  the  wavenumber  wluch  gives  the 
maximum  of  energy. 

In  order  to  take  into  account  wideband  signals  (such  as 
signals  from  local  events),  we  have  used  the  ’Zero-Delay 
Power  Spectrum’  as  suggested  by  Gupta  et  al.  [1989]. 
This  method,  proposed  by  Nawab  et  al.  [1985],  computes 
the  ’Zero-Delay  Covariance  Matrix’  from  the  recorded  sig¬ 
nals,  then  computes  the  wavenumber  spectrum  using  the 
’Maximum  Likelihood  Estimator’.  For  a  plane  wave  with  a 
given  velocity,  the  radial  k-spectrum  in  the  wave  azimuth 
is  the  spectrum  of  the  wave.  The  estimation  of  azimuth  and 
slowness  is  given  by  the  wave  number  which  realizes  the 
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maximum  of  energy.  Another  estimation  of  the  azimuth 
can  be  obtained  by  searching  for  the  azimuth  which  gives 
the  maximum  of  radial  energy  integrated  over  the  radial 
wavenumbers.  In  the  following,  we  refer  to  these  two  kinds 
of  azimuth  estimation  as  the  ‘Maximum  of  Amplitude’  and 
the  'Maximum  of  Radial  Energy’  determinations  respec¬ 
tively. 

The  second  method  is  based  on  the  correlation  of  the 
signals.  In  contrast  to  the  ’k-specmim  Method’  which  com¬ 
putes  the  wave  parameters  by  a  straightforward  process, 
the  ’Correlation  Method’  acts  step  by  step  and  computes 
sequentially  the  arrival-time  differences  for  each  couple  of 
stations,  then  derives  the  wave  parameters  (azimuth  and 
slowness)  using  a  classical  location  method  [e.g.:  Husebye, 
1969].  This  step-by-step  algorithm  allows  us  to  correct  the 
arrival-time  differences  for  plane  wave  model  anomalies 
using  for  example  a  set  of  statistically  determined  station 
corrections. 

The  high  accuracy  measurement  of  the  arrival-time 
differencies  At,(  for  each  couple  of  stations  is  performed  in 
two  steps.  A  rough  estimate  (step  1)  is  yielded  by  the  max¬ 
imum  value  of  the  cross-correlation  function  (maximum 
accuracy  equal  to  the  sampling  interval),  then  a  refined  one 
(step  2)  which  takes  into  account  the  frequency  content  is 
obtained  by  the  slope  of  the  weighted  least-squared  fit  of 
the  cross-spectrum  phase  as  a  function  of  selected  frequen¬ 
cies  which  have  high  signal-to-noise  ratio.  The  coherency 
function  is  used  as  the  weighting  factor.  Nevertheless,  for 
teleseismic  events,  the  slope  of  the  fit  cannot  be  deter¬ 
mined  with  a  good  accuracy  because  of  their  narrow 
frequency-band.  To  remove  this  instability,  we  use  the  fact 
that  the  phase  is  tending  to  zero  with  frequency  and  solve 
the  least-squared  fit  only  for  the  slope. 

For  local  events,  seismic  signals  have  a  higher  fre¬ 
quency  content,  and,  consequently,  in  some  cases  (specially 
for  small  events)  the  cross-correlation  function  has  not  a 
unique  well-defined  maximum.  Two  or  more  delays  might 
give  correlation  factors  close  to  one  another.  This  feature  is 
named  ’ambiguity’. 

To  remove  the'  errors  induced  by  this  ambiguity,  the 
’consistency’  of  the  time-delay  set  is  tested.  For  each  set  of 
stations  i,  j,  k,  let  us  introduce  the  residual  r,jk,  also 
denoted  Chasles  residual: 

rijk  =  +  At*  +  At* 

and  the  error  of  each  time-delay  e,j. 

The  time-delay  set  is  considered  as  ’consistent’  if  the 
triangular  ’Chasles  relationship’: 

f.jt  =  o 

holds  for  all  stations  i,  j,  k.  If  not,  we  have  to  find  the 
e^’s  for  which: 

eij  +  cjk  +  =  rijk 

This  system  is  strongly  underdeteimined:  but  we  can 
translate  it  into  an  integer-number  system  of  equations  by 
assuming  that  all  these  numbers  are  multiples  of  the  sam¬ 
pling  (or  oversampling)  interval.  Furthermore,  we  are  only 
interested  in  the  largest  values  of  the  errors  (i.e.:  large 
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Fig.  2.  Example  of  ambiguity  error  for  a  Pn-wave.  See  text 
for  explanations. 

enough  to  correspond  to  another  maximum  of  the 
corresponding  correlation  function).  Using  linear  program¬ 
ming,  we  can  easily  solve  this  problem,  ana  then  use  the 
largest  values  of  e,j  to  look  for  a  new  consistent  time-deiay 
set. 

Figure  2  shows  an  example  which  describes  clearly  this 
feature.  The  cross-correlation  function  of  a  Pn-wave  is 
shown,  by  dots  for  the  original  digitization  frequency  (50 
Hz),  and  by  a  solid  line  for  an  oversampling  rate  of  8.  The 
triangle  shows  the  maximum  of  this  function,  and  the 
arrow  the  secondary  maximum  which  is  compatible  with 
the  maxima  of  the  other  couples  of  stations  using  Chasles 
relationship.  For  this  event,  table  1  shows  the  results 
(’RMS  Ch’:  root  mean  squared  of  Chasles  residuals.  ’AV 
and  ’AAz’:  velocity  and  azimuth  errors  relative  to  Pn- 
velocity  and  to  true  azimuth  respectively)  in  the  two  cases: 
with  and  without  the  use  of  Chasles  relationship.  It  indi¬ 
cates  clearly  that  valuable  results  are  obtained  when  we  use 
Chasles  relationship  (even  without  oversampling). 

The  Results  for  Teleseismic  Events 

The  two  azimuth  estimation  methods  derived  from  k- 
spectrum  analysis  have  been  applied  to  each  event. 

Figure  3  shows  the  difficulty  of  determinatio"  of  the 
azimuth  using  the  ’Maximum  of  Amplitude’  determination: 
the  solid  curve  represents  the  radial  maximum  of  energy 
versus  azimuth,  which  is  bounded  by  the  k-spectrum  value 
for  k=0.  The  dashed  one  shows  the  variation  of  the  radial 
energy  vertrs  azimuth.  This  last  determination  gives  a 
better  precision  because  of  its  greater  range  of  amplitude. 

The  errors  of  azimuth  determinations  (i.e.:  the 
differences  between  the  k-spectrum  azimuths  and  azimuths 
derived  from  the  USGS  locations)  are  plotted  versus  dis¬ 
tance  on  figure  4.  For  both  kinds  of  determination  (i.e.: 
’Maximum  of  Amplitude’  and  ’Maximum  of  Radial 
Energy’),  the  azimuth  is  not  available  for  distances  greater 
than  80  degrees  (specially  PKP-distances),  but  for  smaller 
distances  the  determination  is  available  and  the  standard 
deviation  is  18  degrees  for  the  ’Maximum  of  Amplitude' 


Table  1 


Use  of  Chales 

Sampling 

RMS  Ch 

AAz 

AV 

Relations  hio 

rate  (Hz) 

(s) 

(°) 

(km/s) 

no 

400 

0.062 

6.2 

1.74 

yes _ 

50 

0.002 

0.0 

0.09 
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Fig.  3.  Variation  of  radial  maximum  amplitude  (solid  line) 
and  radial  energy  (dashed  line)  versus  azimuth. 

determination  and  15  degrees  for  the  ’Maximum  of  Radial 
Energy’  one. 

The  same  dataset  has  also  been  processed  by  the 
’Correlation  Method’.  For  each  event,  the  RMS  value  of 
the  Chasles  residuals  r^  is  less  than  0.003s,  which  gives 
an  evaluation  of  the  arrival-time  differences  measurements. 

Using  tne  rime-delay  sets  of  all  the  event:,  we  have 
localized  them  and  then  computed  station  anomalies  as  th 
mean  values  of  the  residuals.  The  greatest  one  (GM7: 
0.008s)  cannot  be  neglected  compared  with  the  measure¬ 
ment  accuracy:  0.003s. 

Final  results,  taken  into  account  the  station  corrections, 
are  presented  in  figure  5.  The  azimuth  errors  are  plotted 
versus  USGS  azimuth.  It  shows  clearly  a  cosine  depen¬ 
dence  which  might  be  explained  by  a  dipping  structure  of 
the  crust  layers  under  the  mini-array.  After  removal  of  this 
cosine  dependence,  the  standard  deviation  is  reduced  to  9 
degrees. 


*  6=  20° 
♦  6=  17° 


The  Results  for  Regional  Events 

Concerning  the  regional  events,  we  have  only  used  the 
’Correlation  Method’.  The  high  frequency  content  of  the 
signals  (specially  for  the  Pn-phase)  leads  to  wavelengths 
of  the  order  of  the  size  of  the  mini-array.  So.  the  k- 
spectrum  contains  many  relative  maxima  which  have 
amplitudes  very  close  to  each  other. 


Fig.  4.  Azimc'h  errors  (relative  to  USGS  azimuths)  as  a 
function  of  >  ..ince,  using  ’Maximum  of  Radial  Energy’. 
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Method-:  azimuth  error  is  plotted  versus  USGS  azimuth. 
Filled  symbols  are  not  used  for  the  computation  of  cosine 
dependence. 

The  results  obtained  with  the  'Correlation  Method’  are 
summarized  on  figure  6,  which  shows  the  azimuth  errors 
versus  the  LDG  azimuth  taken  as  the  reference.  One  can 
see  again  a  small  cosine  dependance,  which  might  not  be 
significant.  Tire  residual  standard  deviation  is  in  this  case  3 
degrees. 


The  ’Correlation  Method’  used  as  a  Detector 

As  described  above,  Chasles  relationship  is  used  to  test 
the  consistency  of  the  time-delay  set.  When  applied  to 
plane  wave  signals,  the  RMS  of  Chasles  residuals  is 
comparable  to  the  measurement  accuracy  (i.e.:  0.003s). 

If  we  apply  this  process  to  the  noise  recorded  on  the 
stations,  RMS  of  Chasles  residuals  is  about  0.01s  or  more, 
to  be  compared  to  the  0.003s  obtained  for  seismic  signals. 
Consequently,  we  can  use  this  quantity  as  a  discriminant 
between  seismic  signals  and  noise. 

Figure  7  shows  an  example  for  a  local  event.  For 
different  frequency-bands,  the  signals  are  analysed  by 
using  a  moving  rime-window  of  2s-width  .  The  plot  on  the 
top  shows  a  colored  map  of  the  velocity  obtained  for  each 
position  of  the  time-window  (horizontal  axis)  and  of  the 


Pn  phase:  triangles  -  Pg  phase:  crosses 


Fig.  6.  Final  results  for  regional  events  using  ’Correlation 
Method’:  azimuth  error  is  plotted  versus  LDG  azimuth. 
Filled  symbols  are  not  used  for  the  corr  .icn  of  cosine 
dependence. 
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Fig  7.  Example  of  phase  detection  and  identification  by 
using  time-frequency  distribution  of  the  velocity  computed 
only  for  consistent  parts  of  signals  (i.e  Chasles  residuals 
RMS  less  than  0.009s) 

frequency-band  (vertical  axis).  But  the  velocity  color  is 
only  present  if  the  RMS  of  Chasles  residuals  is  less  than 
0009s 

It  clearly  appears  that  only  the  uine-windows 
corresponding  to  the  beginning  of  each  phase  match  this 
criterion  The  influence  of  the  frequency  band  is  also 
clearly  seen:  the  Pn-phase  is  consistent  in  the  high  fre¬ 
quency  domain,  but  Lg  phase  consistency  is  limited  to  low 
frequencies.  Furthermore,  the  velocity  distribution  allows  a 
phase  identification  without  any  ambiguity 

Conclusions 

We  have  tested  two  methods  devoted  to  automatic  loca¬ 
tion  of  both  teleseismic  and  regional  events.  They  are 
based  on  a  major  assumption:  the  very  good  correlation  of 
the  different  records  of  each  event  between  stations.  More¬ 
over,  another  assumption  is  needed  by  the  k-spectrum 
method:  it  assumes  that  a  plane  wave  model  is  valid  for  all 
the  data 

For  teleseisms,  using  the  azimuth  determined  by  the 
'Maximum  of  Radial  Energy’  of  the  k-spectrum,  the 
azimuth  standard  error  is  about  15  degrees  for  distances 
smaller  than  80  degrees.  On  the  opposite,  no  a  priori 
model  is  needed  for  the  correlation  method  This  allows  us 
to  proceed  in  two  steps:  firstly,  a  rough  At,j  evaluation 
(correlation),  and  then  a  refined  one  (phase  spectrum),  with 
the  possibility  of  introducing  station  anomalies  before  loca¬ 
tion.  Nevertheless,  a  careful  study  ’ch  time-delay  set  is 


necessary  to  be  sure  of  its  consistency  (Chasles  relation¬ 
ship).  Using  this  method,  the  standard  enor  on  azimuth  is 
reduced  to  less  than  10  degrees,  after  correction  of  an  evi¬ 
dent  cosine  azimuthal  dependance.  This  dependance  might 
be  due  to  a  dipping  structure  of  the  crustal  layers  under  the 
mini-anay. 

Similarly,  for  regional  events,  the  azimuth  standard 
en-or  is  less  than  3  degrees,  for  a  data  set  including  both 
Pn  and  Pg  phases. 

Finally,  we  have  tested  the  detection  possibilities  of  the 
correlation  method  by  using  a  criterion  based  on  the  con¬ 
sistency  of  the  Chasles  residuals.  It  allows  us  to  detect  and 
to  identify  the  later  arrivals  inside  each  signal,  cuch  as  Sn 
and  Sg  waves  for  regional  events.  Further  studies  will 
investigate  the  potential  of  this  method  in  some  other  cases 
such  as  a  new  implementation  of  the  mini-array  using  10 
stations  and  also  its  application  to  a  larger  network  of 
regional  extension 
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OBJECTIVES 

Among  the  different  methods  used  to  automatically  locate  earthquakes  recorded  on  a  mini-array,  the 
method  based  on  the  broad  band  F-k  analysis  is  the  most  common  (e.g.  :  Mykkeltveit  et  al.,  1990).  However,  this 
method  assumes  that  the  propagation  of  the  considered  wave  front  can  be  modeled  by  a  plane-wave  at  the  scale 
of  the  array.  Furthermore,  in  the  case  of  local  model  variations  (e.g. :  station  anomalies)  the  resulting  precision  is 
very  low  because  of  the  large  wavelength  compared  to  the  array  extension.  To  take  into  account  these  difficulties, 
we  have  tested  another  method  based  on  a  high  precision  determination  of  the  arrival-times,  from  which  the 
event  location  is  derived  by  a  classical  Husebye's  method  which  computes  the  velocity  and  the  azimuth  of  the 
wave.  Each  time-window  is  processed  by  this  method  in  different  frequency  bands  and  activated  only  if  the 
arrival-time  delay  set  is  consistent  Furthermore,  this  consistency  is  used  as  a  signal  detector  and  leads  to  three 
time-frequency  functions,  representing  first  the  consistency,  and  second  the  velocity  and  azimuth  when  they  are 
available.  In  the  case  of  teleseismic  events,  the  event  location  is  strictly  derived  from  the  determinations  of  the 
velocity  and  the  azimuth.  But  for  regional  events  a  phase  identification  is  needed.  This  task  is  performed  by  a 
neural  network  which  uses  the  three  time-frequency  functions  as  inputs  and  which  leads  to  an  estimate  of  each 
phase  occurrence  possibility  as  a  function  of  time. 

Processing  of  a  set  of  28  regional  events  recorded  on  the  5-station  mini-array  worked  by  the 
Laboratoire  de  Detection  et  de  Geophysique  of  the  french  Commissariat  i  lEnergie  Atomique  led  to  the 
following  results : 

-  a  very  precise  determination  of  the  azimuth  :  standard  error  is  less  than  4  degrees; 

-  a  good  phase  identification  capability  (especially  for  P-waves)  which  leads  to  an  estimation  of  the 
distance  with  a  relative  error  lower  than  20%. 


RESEARCH  ACCOMPLISHED 

In  1992,  the  french  Laboratoire  de  Detection  et  de  Gdophysique  of  the  Commissariat  h  TEnergie 
Atomique  has  installed  in  Central  France  a  temporary  mini-array  composed  of  5  vertical  short-period 
seismometers  (Figure  1).  Numerical  signals  are  digitized  at  a  rate  of  50  cycles  per  second  with  a  12-bit  dynamic. 

This  mini-array  has  recorded  more  than  100  teleseismic  events  and  about  28  regional  events  during  its 
4-months  operating  period.  This  experience  has  provided  a  useful  database  for  testing  different  automatic 
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location  methods  (Cansi  et  al.,  1992).  We  have  shown  that  better  results  are  obtained  when  we  use  the  two-step 
correlation-method  which  computes  first  the  arrival-time  differences  with  a  high  precision  (less  than  the 
sampling  interval  0.02s),  and  second  the  azimuth  and  propagation  velocity  corrected  for  statistically  determined 
station  anomalies. 

1)  Data  Processing  : 

For  regional  events,  the  correlation  functions  which  define  the  arrival-time  differences  often  have 
several  relative  maxima  with  comparable  amplitudes.  This  leads  to  an  ambiguity  in  the  arrival-time  computation, 
which  can  be  removed  by  testing  the  consistency  of  the  time-delay  set,  using  the  Chasles  relationship  : 

Aty  +  Atjk  +  Ata  =  0 

Furthermore,  this  consistency  can  be  used  as  a  signal  detector.  When  the  studied  time-window 
contains  a  seismic  signal,  the  Root  Mean  Square  of  the  residuals  of  the  Chasles  relations  is  low  (i.e.  :  less  than 
the  sampling  interval :  it  is  an  estimate  of  the  measurement  accuracy).  On  the  opposite,  when  it  contains  only 
noise,  no  consistency  can  be  found  (i.e.  :  the  RMS  of  Chasles  relations  is  high),  because  of  its  low  correlation  at 
the  scale  of  the  array. 

Then,  for  each  4.5s-time  window  and  for  different  frequency-bands,  we  can  estimate  a  probability  of 
sigral  occurrence,  and,  in  the  case  of  high  probability,  the  corresponding  velocity  and  azimuth.  Some  examples 
of  these  time-frequency  functions  are  displayed  on  figures  3  and  4.  For  each  time  step  and  for  each  frequency- 
band,  the  velocity  is  shown  (grey  scale)  only  when  the  consistency  is  better  than  0.02s. 

We  can  see  clearly  that  most  of  the  seismic  phases  lead  to  consistent  signals  whose  velocity  is  well 
defined  for  all  the  available  frequency-bands.  Nevertheless,  some  cases  are  more  ambiguous  : 

-  a  phase  cannot  be  precisely  recognized  because  the  velocity  is  not  clearly  defined  (see  Sn-phase  on 
figure  4), 

-  a  false  detection  is  obtained  because  a  part  of  the  record  contains  consistent  noise  with  a  velocity 
compatible  with  the  regional  phases  velocity  (see  noise  on  figures  3  and  4  as  an  example). 

Since  these  two  kinds  of  problems  cannot  be  easily  modeled,  we  have  used  a  "learning  system" 
approach  based  on  a  neural  network  to  identify  the  different  phases  of  each  event  without  ambiguity. 

2)  Phase  Identification  : 

For  neuromimetic  applications,  the  classical  programming  efforts  are  transposed  to  the  determination 
of  the  various  authorised  degrees  of  freedom  described  as  follows  : 

-  the  data  structure  :  the  first  step  is  to  extract  from  the  database  the  information  that  are  strictly 
necessary  for  phase  identification.  Furthermore,  those  data  have  to  be  invariant  by  translation,  rotation  and 
dilatation,  which  precludes  the  analysis  on  a  variable  period.  The  solution  we  chose  is  thus  to  present  as  inputs 
and  for  each  time  sample  the  signal  velocity  for  5  frequency  bands  (  from  0  to  12  Hz).  The  data  whose 
corresponding  Chasles  RMS  is  greater  than  an  arbitrary  threshold  (l.e. :  0.02  s)  are  set  to  0. 
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-  the  network  structure  :  we  only  used  multi-layered  perceptrons,  with  a  sigmoid  tranfer  function.  They 
are  indeed  capable  of  building  complex  decision  hyper-volumes  in  the  hyper-space  of  the  inp-'.t  data,  thus 
realising  an  accurate  classifier.  Several  tests  led  to  choose  a  2  hidden-layer  perception.  The  complexity  of  the 
system  is  due  to  the  high  non-linearity  of  the  problem. 

-  the  learning  function  :  we  chose  as  a  learning  function  the"back  propagation  with  momentum" 
method,  which  uses  a  gradient  method  to  minimize  the  quadratic  error  between  the  expected  and  the  observed 
results.  Despite  a  long  computing  time,  this  allows  a  reliable  and  accurate  learning  convergence. 

-  the  example  database  :  it  was  made  of  all  the  available  events,  excepted  3  of  them  on  which  the  tests 
were  made.  In  order  to  avoid  incoherencies,  the  arrival  time  of  each  phase  was  picked  on  the  time-frequency 
diagrams  ;  a  phase  is  thus  declared  present  over  the  whole  time-window  between  its  arrival  time  and  the  arrival 
time  of  the  following  phase.  Each  sample  is  presented  50000  times  in  a  random  order. 

-  the  network  topology  :  the  number  of  nodes  in  the  hidden  layers  was  determinated  empirically. 
Several  networks  were  designed :  in  order  to  avoid  over-training,  we  chose  for  each  phase  the  simplest  one  which 
did  not  degrade  the  results.  Finally,  the  Pn  phase  requires  24+9  units,  the  Pg  phase  20+6,  the  Sn  phase  16+6  and 
the  Sn  phase  20+6  (  Figure  2  ). 

All  the  designing  learning  and  testing  operations  were  performed  using  the  neural  simulator  SNNS, 
developed  by  the  Stuttgart  University.  The  middle  diagrams  of  figures  3  and  4  show  the  neural  outputs  as  a 
function  of  time  for  two  events  which  were  not  in  the  learning  database. 

3)  Distance  estimation  : 

In  order  to  remove  the  last  false  detections  due  to  consistent  noise,  a  post  data  processing  is  needed  to 
test  the  consistency  of  the  results  on  the  whole  signal,  by  using  the  azimutal  information  as  describe  as  follows  : 

-  the  first  step  is  to  compute  an  azimuth  histogram  with  the  possibilitiy  fonctions  previously 
determined  and  to  to  choose  the  more  probable  20°  wide  interval.  The  average  azimuth  can  then  be  calculated. 

-  the  second  step  is  to  refine  this  approximation  by  determining  the  most  probable  5°  wide  sub-interval 
for  each  phase.  For  each  time  sample,  the  probability  of  existence  of  a  phase  is  set  to  0  if  the  corresponding 
azimuth  is  out  of  those  sub-intervals. 

The  bottom  diagrams  in  figures  3  and  4  show  the  phase  characterization  curves  after  this  azimuthal 

filtering. 

At  this  step,  in  the  function  describing  the  possibility  of  occurrence  of  each  phase  as  a  function  of 
time,  all  the  information  which  do  not  belong  to  the  detected  event  have  been  removed.  The  last  step  -  the 
distance  evaluation  -  can  now  be  performed. 

Each  function  is  differentiated  to  allow  an  easy  detection  of  each  phase  by  identifying  the  times  where 
the  derivative  is  greater  than  an  arbitrary  threshold  (i.e. :  0.7).  When  only  two  phases  are  detected,  the  resulting 
distance  is  computed  using  the  two  times  for  which  the  product  of  the  corresponding  possibility  functions  is  at  its 
maximum  When  more  than  two  phases  are  detected,  a  least-squares  estimation  of  the  distance  is  performed  for 
each  set  of  possible  arrival-times.  The  best  one  is  retained  as  the  event  distance. 
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Figure  1  :  Network  map.  It  is  centred  on  the  point  45.74 11N-2.0133E. 
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Figure  2  :  Desien  of  the  neural  networks. 
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Figure  3  :  Example  of  an  event  at  the  different  processing  levels  :  the  time-frequency  plot  of  the  velocity  (top), 
the  4  phase  probability  functions  (middle j, ^and  the  same  after  azimuth  filtering  (bottom). 
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OBJECTIVE 


The  French  LDG  (Laboratoire  de  Geophysique)  has  been  installing  and 
operating  an  experimental  mini-array  in  Central  France  (GIAT)  for  seismic 
Verification  and  Monitoring  purposes,  since  March,  1993.  The  objective  of  the 
work  presented  here  was  to  adapt  and  use  the  software  developed  at  Norsar,  in 
order  to  assess  the  detection  capabilities  of  the  GIAT  array,  both  for  regional  and 
teleseismic  events. 


RESEARCH  ACCOMPLISHED 


Data  acquisition  and  processing: 

The  French  GIAT  mini-array  presently  consists  of  9  vertical  seismometers, 
and  one  3-components  station.  It's  size  is  about  3  kilometres  (figure  1).  The  data 
are  digitised  at  a  rate  of  50  samples  a  second,  and  are  stored  locally  on  a  personal 
computer  in  files  containing  each  about  17  mn  of  data.  Once  a  file  is  completed,  it 
is  sent  to  the  distant  processing  centre  via  phone,  where  it  is  automatically 
written  to  a  Norsar  type  diskloop.  The  detection  software  recognises  the  arrival  of 
new  data,  computes  a  set  of  filtered  beams,  and  then  uses  a  classical  STA/LTA 
detector  on  each  beam. 

Detection  configuration: 

The  results  reported  here  were  obtained  over  a  57  day  period,  between 
May  5  and  June  30,  1993.  During  this  period,  a  set  of  83  beams  was  deployed,  at 
different  velocities  (for  regional  phases)  and  for  many  azimuths.  These  beams 
were  also  attached  to  different  sub-arrays,  which  were  empirically  and 
qualitatively  determined  according  to  the  frequency  content  and  correlation 
characteristics  of  the  different  wave^pes  and  the  background  noise  (Kvaerna, 


1989).  STA/LTA  thresholds  were  between  4.4  and  7.1  for  coherent  beams,  and 
between  3.6  and  4.4  for  incoherent  beams.  No  horizontal  beams  could  be  used, 
since  only  one  3  components  station  was  available. 

Detection  results: 

During  the  period  considered,  an  average  of  161  detections  was  observed 
each  day.  Deviations  from  this  mean  can  be  very  high  (figure  2),  due  to  transient 
noise  occurring  on  one  or  more  stations  (it's  origin  can  be  cultural, 
meteorological,  traffic,  etc.).  The  smallest  numbers  of  detections  are  sometimes 
observed  on  Saturdays  or  Sundays,  but  this  is  not  always  true. 

As  a  reference  for  estimating  detection  capabilities,  we  used  the  LDG 
regional  and  teleseismic  bulletins,  which  are  obtained  from  the  whole  LDG 
network  (about  40  stations  covering  France),  and  which  are  available  and 
distributed  once  a  week. 

Only  events  for  which  both  a  location  and  a  magnitude  were  reported  in 
these  bulletins  were  used  for  the  comparisons.  Thus,  184  regional  events  and  186 
teleseismic  events  were  selected,  i.e.  an  average  of  about  3  regional  events  and  3 
teleseismic  events  every  day. 

An  event  was  declared  to  be  detected  by  the  mini-array  if  at  least  one 
phase  was  detected  by  at  least  one  beam.  Note  that  this  does  NOT  mean  that 
such  an  event  could  be  located  or  identified  automatically  by  the  mini-array. 

Results  indicate  that  less  than  5  %  of  all  detections  can  be  associated  with 
natural  seismic  events.  Another  5  to  10  %  can  be  attributed  to  artificial  sources 
such  as  quarry  or  rock  blasts.  In  fact,  from  Mondays  to  Fridays,  the  mini-array 
detected  an  average  of  5  quarry-blasts  per  day  during  working  hours.  Thus, 
nearly  90  %  of  all  detections  were  never  attributed  to  a  known  seismic  source. 

Detectability  of  regional  events: 

Figures  3  and  4  display  the  preliminary  results  obtained  for  regional 
events.  It  can  be  concluded  from  these  figures  that  the  75%  confidence  threshold 
for  detection  is  about  Ml  =  3.4  (Ml  is  the  LDG  local  magnitude). 

This  value  is  about  0.7  magnitude  units  higher  than  the  value  obtained  for 
the  whole  LDG  network,  which  can  be  estimated  as  Ml  =  2.7  (LDG  local 
magnitude),  for  distances  up  to  6  to  8  degrees. 

It  is  instructive  to  compare  this  threshold  with  those  reported  for  the 
fennoscandian  arrays  (Mykkeltveit  et  al.  1990;  Uski,  1990),  which  are  close  to  Ml 
=  2.3  to  2.7  (where  Ml  is  the  local  magnitude  computed  for  Norway)  at  the  90% 
confidence  level. 

Thus,  the  difference  between  the  detection  capabilities  at  both  sites  is 
presently  close  to  one  unit  of  local  magnitude.  This  is  apparently  due  to  the 
different  calibrations  used  for  these  local  magnitudes,  which  should  be  further 
investigated  in  the  future.  Also,  there  is  about  a  factor  two  difference  in  Lg  wave 
attenuation  coefficients  between  Norway  and  France  (Alsaker  et  al.  1991;  Plantet 
et  al.  1991). 


Detectability  for  teleseismic  events: 
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Figures  5  and  6  display  the  results  obtained  for  the  detectability  of 
teleseismic  events.  The  75%  confidence  threshold  can  presently  be  estimated  at 
mb  =  4.9,  which  is  close  to  the  whole  LDG  network  threshold,  though  the  LDG 
bulletin  can  probably  no  longer  be  used  as  a  reference  in  the  teleseismic  case. 
Anyway,  only  one  event  with  a  magnitude  above  5.4  was  missed. 

These  results  are  close  to  those  reported  for  GERESS  (Gestermann  et  al. 

1991). 

Regional  phase  identification : 

Regional  phases  observed  in  central  France  are  mainly  Pn,  Pg,  Sg  and  Lg, 
while  Sn  is  rarely  detected.  For  Fennoscandian  arrays,  Norsar  uses  both  phase 
velocity  (obtained  from  the  f-k  analysis)  and  polarisation  results  in  order  to 
identify  regional  phases.  However,  this  is  more  difficult  for  GLAT,  since  only  one 
(very  noisy)  3  components  station  is  available.  So,  presently,  we  only  use 
velocities  for  identifying  regional  phases. 

An  analysis  of  38  such  phases,  carefully  checked  by  the  analyst,  led  to  the 
following  preliminary  rules  (velocities  are  given  in  km/sec): 

if  VEL  <3.1  then  phase  is  Noise 
else:  if  VEL  <  5.8  then  phase  is  S 

if  VEL  >3.1  and  VEL  <  4.2  then  phase  is  Lg 
if  VEL  >  4.2  and  VEL  <  5.8  then  phase  is  Sn 
else:  if  VEL  >  5.8  then  phase  is  P 

if  VEL  >  5.8  and  VEL  <  7.4  then  phase  is  Pg 
if  VEL  >  7.4  and  VEL  <  10.5  then  phase  is  Pn 
else:  if  VEL  >  10.5  then  phase  is  teleseismic 

The  main  ambiguity  which  may  appear  is  for  Sg-Lg  versus  Sn 
discrimination,  for  velocities  between  4.2  and  5.2.  More  experience  is  now  needed 
before  evaluating  the  performances  of  these  phase  identification  rules. 


CONCLUSIONS  AND  RECOMMENDATIONS 


Despite  the  fact  that  the  GIAT  mini-array  consists  of  only  10  stations, 
detection  results  for  both  regional  and  teleseismic  events  are  encouraging:  the  75 
%  confidence  threshold  is  ML=3.4  for  regional  and  mb=4.9  for  teleseismic  events. 

Such  figures  could  only  be  obtained  due  to  the  high  false  detection  rate 
(nearly  90%)  that  was  allowed.  It  is  expected  that  these  false  detections  will  be 
eliminated  during  the  phase  association  procedure  included  in  the  RONAPP 
software,  and  that  no  or  little  spurious  events  will  be  found  in  the  final  automatic 
bulletin. 

We  are  presently  trying  to  further  reduce  the  STA/LTA  thresholds,  to 
ameliorate  the  spike  elimination  algorithm,  and  to  increase  the  total  number  of 
beams.  We  plan  to  run  with  an  average  of  perhaps  400  detections  a  day,  among 
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which  less  than  5%  would  be  associated  with  real  seismic  events.  We  hope  that 
this  will  further  reduce  the  detection  thresholds  by  about  0.2  to  0.3  magnitude 
units.  However,  it  is  not  yet  sure  weather  regional  events  with  such  low 
magnitudes  (say  3.1  local  LDG  magnitude)  could  be  automatically  located,  since 
for  that  purpose  at  least  two  phases  must  be  detected. 

A  longer  period  of  observation,  and  comparisons  with  other  reference 
bulletins  (for  teleseismic  events)  should  confirm  these  results. 

Of  course,  adding  more  stations,  especially  3-component  stations,  would 
also  improve  the  results. 

The  next  step  is  now  to  investigate  in  more  details  the  phase  identification 
and  location  capabilities  of  the  GIAT  mini-array.  This  will  be  done  using  both  the 
classical  RONAPP  algorithm  provided  by  Norsar,  and  a  novel  approach  for  phase 
picking  and  classification  based  on  coherency  analysis  and  neural  networks 
(Cansi  et  al.,  1993). 

Finally,  before  distributing  automatic  bulletins,  it  will  be  necessary  to 
address  the  important  issue  of  discriminating  between  quarry  blasts  and  natural 
seismic  events. 
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FIGURE  1.  Location  and  geometry  of  the  French  GIAT  mini-array  in  central 
France.  The  array  consists  of  9  vertical  seismometers,  and  one  "central"  3- 
components  station  (GM5).  GM5  is  located  at  45°44  N  and  2°01  E. 


FIGURE  2.  Number  of  detections  per  day  at  GIAT,  from  May  5  to  June  30,  1993, 
deploying  a  set  of  83  beams.  Mean  is  161.  Deviations  from  the  mean  by  more 
than  a  factor  2  are  not  rare.  The  effect  (jf^eekends  does  not  appear  clearly. 
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FIGURE  3.  Regional  even,  detectability  at  GIAT.  The  Y-axis  displays  the 
cumulated  number  of  detected  events,  whose  local  LDG  magnitude  is  higher  than 
the  value  specified  on  the  X-axis.  Reference  is  the  LDG  regional  bulletin. 
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FIGURE  4.  Regional  event  detectability  at  GIAT.  Reference  is  the  LDG  regional 
bulletin.  The  75%  confidence  threshold  for  detection  is  3.4.  Up  to  a  distance  of  8 
degrees,  no  event  with  local  magnitude  greater  than  3.2  was  missed. 
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FIGURE  5.  Telfeseismic  event  detectability  at  GIAT.  The  Y-axis  displays  the 
cumulated  number  of  detected  events,  whose  body  wave  magnitude  mb  is  higher 
than  the  value  specified  on  the  X-axis.  Reference  is  the  LDG  teleseismic  bulletin. 
Theoretical  relation  is  linear  with  a  slope  equal  to  -1.  Deviation  from  linearity  is 
observed  for  an  mb  of  about  4.8  and  points  out  the  detection  threshold. 
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FIGURE  6.  Teleseismic  event  detectability  at  GIAT.  Reference  is  the  LDG 
teleseismic  bulletin.  The  75%  confidence  threshold  for  detection  is  4.9.  Up  to  a 
distance  of  60  degrees,  no  event  with  pjtj?  greater  than  4.9  was  missed. 
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